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EXECUTIVE  SUMMARY 

The  passive  location  of  both  hostile  and  friendly  electromagnetic  emitters  has  been 
an  important  capability  for  the  war-fighter,  for  law  enforcement,  and  in  search  and  rescue 
operations.  The  characteristics  of  an  emitted  signal  once  received  by  several  sensors  in  an 
operating  environment  can  be  exploited  to  passively  locate  an  emitter.  Two  of  the 
characteristics  that  can  be  measured  by  separate  receivers  of  such  a  signal  are:  (1)  the 
differences  in  Doppler  shifts  (if  any)  in  the  signal's  frequency  and  (2)  the  difference  in  the 
time  that  the  signal  arrived  at  the  sensors.  These  measurements  are  known  as  Frequency 
Difference  of  Arrival  (FDOA)  and  Time  Difference  of  Arrival  (TDOA).  In  the  presence  of 
moving  sensors,  the  sensors'  receivers  measure  the  frequency  of  arrival  (FOA)  of  the 
signal.  The  difference  of  two  FOA's  provides  an  FDOA  which  yields  a  surface  called  an 
isodop,  in  this  case  meaning  equally  differenced  Doppler  data,  which  contains  the  locus  of 
all  points  on  which  the  emitter  could  lie  given  the  Doppler  information.  Similarly,  the 
difference  of  a  single  receiver's  time  of  arrival  (TO  A)  with  another  TO  A  yields  a  TDOA 
which  describes  a  surface  known  as  an  isochron,  in  this  case  meaning  equally  differenced 
time,  which  contains  the  locus  of  all  points  on  which  the  emitter  could  lie  given  the  TDOA 
information. 

In  the  past,  the  measurement  of  such  quantities  in  the  tactical  environment  has 
been  difficult  due  to  their  sensitivities  to  timing  errors  between  any  two  sensors.  Further, 
when  trying  to  use  these  measurements  for  the  geolocation  of  an  emitter,  errors  in  the 
sensor's  own  location  and  velocity  measurements  compound  with  the  signal  processing 
errors  to  provide  unreliable  emitter  fixes.  Historically,  the  Navy  has  relied  upon  long 
range,  shore-based  Fligh  Frequency  Direction  Finding  (HFDF)  networks.  These  networks 
use  large  antennas  to  measure  the  angle  of  arrival  of  a  signal  and  cross  these  bearings  over 
long  distances  to  triangularize  the  emitter's  location.  Given  the  recent  losses  of  Navy 
overseas  assets  and  that  more  and  more  emitters  in  the  tactical  environment  operate  in  the 
Very  High  Frequency  (VHF)  and  Ultra  Fligh  Frequency  (UHF)  range,  the  ability  to 
perform  emitter  geolocation  is  diminishing.  As  a  result,  the  Naval  Security  Group  Support 
Activity  (NSGSA)  contracted  with  the  Applied  Research  Lab  at  the  University  of  Texas  at 
Austin  (ARL:UT)  to  develop  an  affordable,  low-risk  TDOA  geolocation  system  using 
commercial  off  the  shelf  (COTS)  and  government  off  the  shelf  (GOTS)  technologies. 
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ARL:UT  responded  by  the  development  of  the  Carry-On  Multi-Platform  Global 
Positioning  System  Assisted  TDOA  System. 

When  using  TDOA  techniques,  errors  in  the  emitter  geolocation  problem  is 
essentially  related  to  three  components:  the  location  of  the  observer  relative  to  the  emitter 
(the  "geometry"),  signal  timing,  and  observer  position  measurement.  The  emergence  of 
Global  Positioning  System  (GPS)  technology  offers  capabilities  which  can  greatly  enhance 
TDOA  applications  by  reducing  measurement  error  components.  Specifically,  GPS  allows 
for  dramatic  reductions  in  the  observer  position  error  and  sample  timing  between 
observers.  In  addition,  the  advantages  of  using  such  commercially  available  and 
government  "off-the-shelf'  technology  allows  the  system  to  be  placed  on  a  variety  of 
platforms  and,  thus,  could  make  the  GPS-Assisted  TDOA  Geolocation  System  a  highly 
desirable  Multi-Platform  tool  in  the  "Joint  Warfare"  environment. 

The  measurement  of  the  FDOA  is  done  simultaneously  when  measuring  the  TDOA 
in  the  formulation  used  by  the  ARL:UT  Test  System.  While  the  initial  test  and 
development  of  the  ARL:UT  system  focused  on  the  use  of  only  TDOA  measurements  for 
geolocation,  the  combined  use  of  TDOA  and  FDOA  measurements  for  geolocation  has 
become  possible  due  to  the  GPS  technology  in  the  ARL:UT  system.  By  providing  more 
information  to  the  geolocation  solution,  combined  TDOA  and  FDOA  measurements  can 
expand  the  size  of  the  state  of  the  observer — providing  velocity  information  as  well  as  the 
emitter's  location.  In  addition  to  providing  velocity  information,  it  has  been  postulated 
that  when  poor  geometry  exists  between  observers,  FDOA  can  more  tightly  resolve  the 
solution  of  the  geolocation  fix  than  TDOA  alone.  Further,  when  the  emitter  is  believed  to 
have  zero  velocity,  FDOA  can  provide  the  emitter's  three  dimensional  location  using  only 
three  observers  instead  of  four  by  fixing  the  solution's  velocity  state  to  zero. 

Regardless  of  which  method  or  algorithm  is  used,  the  measurement  of  the  TDOA 
and  the  FDOA  is  a  statistical  process,  built  upon  standard  estimation  assumptions.  As 
such,  any  measurement  of  either  component  can  be  viewed  as  a  random  variable  related  to 
the  signal,  the  environment,  and  the  accuracy  and  precision  of  the  whole  system  that 
performs  the  measurement.  Through  the  use  of  a  Java-based,  Stochastic  TDOA  and 
FDOA  Simulation  (JBSTAFSim),  this  thesis  explores  the  impact  on  the  geolocation  of  an 
emitter  in  the  presence  of  such  measurement  errors.  Given  an  estimate  of  the  distribution 
of  these  measurement  errors,  JBSTAFSim  can  simulate  the  impact  on  a  TDOA  or  mixed 
TDOA/FDOA  geolocation  solution  within  a  given  "Joint  Warfare"  scenario  of  observers. 
In  particular,  JBSTAFSim  allows  the  user  to:    adjust  the  magnitude  and  correlation  of 
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measurement  error  variance;  adjust  the  believed  measurement  error  variance  used  by  the 
geolocator;  dictate  the  network  geometry  and  types  of  observers;  and  elect  to  choose  a 
batch  processed  geolocator  which  incorporates  previous  geolocation  information  into  the 
current  solution  to  solve  for  the  emitter's  state. 

Through  simulation,  this  thesis  shows  that  the  mixed  T/FDOA  solution 
stochastically  dominates  the  TDOA  solution  in  two-dimensional  fix  scenarios.  While  no 
stochastic  dominance  of  either  solution  is  shown  in  the  three-dimensional  fix  case,  the 
author  shows  that  the  accuracy  of  the  three-dimensional  problem  is  related  to  the  sensor 
network  geometry.  Further,  if  a  three-dimensional  fix  is  required,  it  is  shown  that  robust 
sensors  like  satellites  should  be  used  to  improve  the  sensor  network  to  target  geometry. 
Finally,  since  the  location  of  the  target  relative  to  the  sensors  is  generally  not  known,  the 
author  demonstrates  the  need  to  define  a  method  to  rate  the  sensor  network  geometry 
given  possible  target  locations. 

These  results  can  be  viewed  from  three  perspectives.  The  first  view  is  from  the 
perspective  of  the  designer  of  the  system  who  wishes  to  understand  the  sensitivities  of 
measurement  errors  and  wishes  to  improve  upon  these  measurements  by  the  system  or 
take  them  into  account  in  the  geolocation  process.  The  next  view  is  from  the  warrior,  or 
user  of  the  system,  who  would  be  employing  this  system  to  localize  a  hostile  or  friendly 
emitter.  The  warrior  is  concerned  with  the  quality  of  the  information  provided  by  the 
system.  As  such,  the  warrior  needs  to  know  how  well  the  system  can  perform  and  what 
he  or  she  can  do  with  respect  to  the  sensor  network  to  improve  the  fix  of  the  target. 
Finally,  the  last  perspective  is  from  that  of  the  activities  which  fund  the  designer  to  create 
such  systems  for  the  warrior.  This  activity  is  concerned  with  system  performance  for  the 
amount  invested.  Given  several  functional  areas  that  system  designers  could  pursue  to 
develop  or  improve,  the  simulation  can  rate  the  impact  of  this  additional  or  improved 
information  on  the  system's  performance  and  describe  how  this  impacts  the  warfare 
commander's  utility  for  such  information. 
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I.  INTRODUCTION 

A.  BACKGROUND 

The  passive  location  of  both  hostile  and  friendly  electromagnetic  emitters  has  been 
an  important  capability  for  the  war-fighter,  for  law  enforcement,  and  in  search  and  rescue 
operations.  The  characteristics  of  an  emitted  signal  once  received  by  several  sensors  in  an 
operating  environment  can  be  exploited  to  passively  locate  an  emitter.  Two  of  the 
characteristics  that  can  be  measured  by  separate  receivers  of  such  a  signal  are:  ( 1 )  the 
differences  in  Doppler  shifts  (if  any)  in  the  signal's  frequency  and  (2)  the  difference  in  the 
time  that  the  signal  arrived  at  the  sensors.  These  measurements  are  known  as  Frequency 
Difference  of  Arrival  (FDOA)  and  Time  Difference  of  Arrival  (TDOA).  In  the  presence  of 
moving  sensors,  the  sensors'  receivers  measure  the  frequency  of  arrival  (FOA)  of  the 
signal.  The  difference  of  two  FOA's  provides  an  FDOA  which  yields  a  surface  called  an 
isodop,  in  this  case  meaning  equally  differenced  Doppler  data,  which  contains  the  locus  of 
all  points  on  which  the  emitter  could  lie  given  the  Doppler  information  (see  Figure  1-1). 
Similarly,  the  difference  of  a  single  receiver's  time  of  arrival  (TOA)  with  another  TOA 
yields  a  TDOA  which  describes  a  surface  known  as  an  isochron,  in  this  case  meaning 
equally  differenced  time,  which  contains  the  locus  of  all  points  on  which  the  emitter  could 
lie  given  the  TDOA  information  (see  Figure  1-2). 


observers 


Figure  1-1  FDOA  Isodops  (From  Ref  1) 


Figure  1-2  TDOA's  of  Three  Observers  (From  Ref.  3) 


In  the  past,  the  measurement  of  such  quantities  in  the  tactical  environment  has 
been  difficult  due  to  their  sensitivities  to  timing  errors  between  any  two  sensors.  Further, 
when  trying  to  use  these  measurements  for  the  geolocation  of  an  emitter,  errors  in  the 
sensor's  own  location  and  velocity  measurements  compound  with  the  signal  processing 
errors  to  provide  unreliable  emitter  fixes.  Historically,  the  Navy  has  relied  upon  long 
range,  shore-based  High  Frequency  Direction  Finding  (HFDF)  networks.  These  networks 
use  large  antennas  to  measure  the  angle  of  arrival  of  a  signal  and  cross  these  bearings  over 
long  distances  to  triangularize  the  emitter's  location.  Given  the  recent  losses  of  Navy 
overseas  assets  and  that  more  and  more  emitters  in  the  tactical  environment  operate  in  the 
Very  High  Frequency  (VF£F)  and  Ultra  High  Frequency  (UHF)  range,  the  ability  to 
perform  emitter  geolocation  is  diminishing.  As  a  result,  "the  Naval  Security  Group 
Support  Activity  (NSGSA)  contracted  with  the  Applied  Research  Lab  at  the  University  of 
Texas  at  Austin  (ARLUT)  to  develop  an  affordable,  low- risk  TDOA  geolocation  system 
using  commercial  off  the  shelf  (COTS)  and  government  off  the  shelf  (GOTS) 
technologies"[Ref  1].  ARL:UT  responded  by  the  development  of  the  Carry-On  Multi- 
Platform  Global  Positioning  System  Assisted  TDOA  System. 

When  using  TDOA  techniques,  errors  in  the  emitter  geolocation  problem  is 
essentially  related  to  three  components:  the  location  of  the  observer  relative  to  the  emitter 
(the  '"geometry"),  signal  timing,  and  observer  position  measurement.  The  emergence  of 
Global  Positioning  System  (GPS)  technology  offers  capabilities  which  can  greatly  enhance 
TDOA  applications  by  reducing  measurement  error  components.  Specifically,  GPS  allows 
for  dramatic  reductions   in   the   observer  position   error   and   sample   timing  between 


observers.  In  addition,  the  advantages  of  using  such  commercially  available  and 
government  "off-the-shelf  technology  allows  the  system  to  be  placed  on  a  variety  of 
platforms  and,  thus,  could  make  the  GPS-Assisted  TDOA  Geolocation  System  a  highly 
desirable  Multi-Platform  tool  in  the  "Joint  Warfare"  environment  as  illustrated  in  Figure 
1-3. 


Figure  1-3  GPS  Assisted  Geolocation  System  in  Joint  Warfare  (From  Ref  2) 


Through  the  use  of  GPS  technology  and  differential  techniques  combined  with  Rubidium 
standard  oscillators,  the  total  theoretical  signal  timing  error  is  approximately  25  nano- 
seconds which,  when  multiplied  by  the  speed  of  light,  converts  to  7.5  meters  in 
distance.  [Ref.  1] 

The  measurement  of  the  FDO A  is  done  simultaneously  when  measuring  the  TDOA 
in  the  formulation  used  by  the  ARL:UT  Test  System.  While  the  initial  test  and 
development  of  the  ARL:UT  system  focused  on  the  use  of  only  TDOA  measurements  for 
geolocation,  the  combined  use  of  TDOA  and  FDOA  measurements  for  geolocation  has 
become  possible  due  to  the  GPS  technology  in  the  ARL:UT  system.  By  providing  more 
information  to  the  geolocation  solution,  combined  TDOA  and  FDOA  measurements  can 
expand  the  size  of  the  state  of  the  observer — providing  velocity  information  as  well  as  the 
emitter's  location.  In  addition  to  providing  velocity  information,  it  has  been  postulated 
that  when  poor  geometry  exists  between  observers,  FDOA  can  more  tightly  resolve  the 
solution  of  the  geolocation  fix  than  TDOA  alone.   Further,  when  the  emitter  is  believed  to 


have  zero  velocity,  FDOA  can  provide  the  emitter's  three  dimensional  location  using  only 
three  observers  instead  of  four  by  fixing  the  solution's  velocity  state  to  zero. 


B.  PURPOSE 

Regardless  what  method  or  algorithm  used,  the  measurement  of  the  TDOA  and  the 
FDOA  is  a  statistical  process,  built  upon  standard  estimation  assumptions.  As  such,  any 
measurement  of  either  component  can  be  viewed  as  a  random  variable  related  to  the 
signal,  the  signal's  environment,  and  the  accuracy  and  precision  of  the  whole  system  that 
performs  the  measurement.  Through  the  use  of  a  Java-based,  Stochastic  TDOA  and 
FDOA  Simulation  (JBSTAFSim),  the  author  explores  the  impact  on  the  geolocation  of  an 
emitter  in  the  presence  of  such  measurement  noise.  Given  an  estimate  of  the  variance  of 
these  measurement  errors,  JBSTAFSim  can  simulate  the  impact  on  a  TDOA  or  mixed 
TDOA/FDOA  geolocation  solution  within  a  given  "Joint  Warfare"  scenario  of  observers. 
In  particular,  JBSTAFSim  allows  the  user  to:  adjust  the  magnitude  and  correlation  of 
measurement  error  variance,  adjust  the  believed  measurement  error  variance  used  by  the 
geolocator;  dictate  the  network  geometry  and  types  of  observers;  and  elect  to  choose  a 
batch  processed  geolocator  which  incorporates  previous  geolocation  information  into  the 
current  solution  to  solve  for  the  emitter's  state. 

These  results  can  be  viewed  from  three  perspectives.  The  first  view  is  from  the 
perspective  of  the  designer  of  the  system  and  wishes  to  understand  the  sensitivities  of 
measurement  errors  and  wishes  to  improve  upon  these  measurements  by  the  system  or 
take  them  into  account  in  the  geolocation  process.  The  next  view  is  from  the  warrior,  or 
user  of  the  system,  who  would  be  employing  this  system  to  localize  a  hostile  or  friendly 
emitter.  The  warrior  is  concerned  with  the  quality  of  the  information  provided  by  the 
system.  As  such,  the  warrior  needs  to  know  how  well  the  system  can  perform  and  what 
he  or  she  can  do  with  respect  to  the  sensor  network  to  improve  the  fix  of  the  target. 
Finally,  the  last  perspective  is  from  that  of  the  activities  which  fund  the  designer  to  create 
such  systems  for  the  warrior.  This  activity  is  concerned  with  system  performance  for  the 
amount  invested.  Given  several  functional  areas  that  system  designers  could  pursue  to 
develop  or  improve,  the  simulation  can  rate  the  impact  of  this  additional  or  improved 
information  on  the  system's  performance  and  describe  how  this  impacts  the  warfare 
commander's  utility  for  such  information. 


C.  ORGANIZATION 

The  next  chapters  discuss  briefly:  the  configuration  of  the  ARL:UT  test  system, 
how  the  TDOA  and  FDOA  measurements  are  made,  and  how  these  measurements  enter 
into  the  geolocation  process.  Some  of  the  detailed  theory  and  methods  that  JBSTAFSim 
uses  can  be  found  in  Appendixes  A-C.  Chapters  V  and  VI  provide  an  overview  of 
JBSTAFSim  and  a  description  of  the  simulated  environment.  Chapter  VII  presents  the 
results  of  the  simulation  of  TDOA  and  mixed  TDOA/FDOA  geolocation  methods. 
Chapter  VIII  offers  recommendations  based  on  the  results  presented  in  the  previous 
chapter. 


H.  ARL:UT  SYSTEM  OVERVIEW 

The  ARL:UT  system  is  capable  of  receiving  a  frequency  in  the  HF,  VHF,  or  UHF 
frequency  bands  Only  one  voice  grade  channel  compatible  with  existing  Navy  hardware 
for  communication  between  sensors  is  required.  Further,  the  system  is  designed  to  be 
compatible  with  the  Unified  Build  (UB)  environment  and  interface  with  and  display  its 
results  graphically  on  the  Joint  Maritime  Command  and  Information  System  (JMCIS). 
Each  observer  updates  its  location  with  a  GPS  receiver.  The  receivers  have  a  one  pulse- 
per-second  (pps)  output  which  is  used  to  trigger  sampling  of  the  incoming  signal  of 
interest.  The  accuracy  of  the  GPS  clocks  drives  each  receivers'  output  of  the  one  pps 
pulse  to  be  within  20  nano  seconds  of  each  other.  The  sampler  at  each  observer,  then,  is 
phase-locked  to  the  one  pps  pulse  to  within  five  nano  seconds  making  the  samples 
between  observers  to  be  within  25  nano  seconds  [Ref.  20].  As  depicted  in  Figure  2-1,  the 
system  network  consists  of  one  master  sensor  and  a  few  slave  sensors.  In  the  evenly 


Figure  2-1  System  Overview  (From  Ref.  20) 


determined  case,  three  sensors  are  required  for  a  two-dimensional  geolocation  and  four  for 
a  three-dimensional  fix.  Once  given  the  frequency  of  interest  from  traditional  detection 
means,  the  master  tunes  its  receiver  to  that  frequency  and  orders  the  slaves  to  tune  to  the 
same  frequency.  Each  observer  stores  the  incoming  signal  as  buffered  data.  Since  only 
one  communications  channel  is  required  to  be  available,  the  ARL.UT  system  uses  a 
distributed  processing  system  to  make  the  most  efficient  use  of  the  single  communications 
link.  The  master  determines  a  400  milli  second  portion  that  contains  the  most  energy  and 
computes  the  Fast  Fourier  Transform  (FFT)  of  this  sample.   This  FFT  is  then  sent  to  the 


slaves  with  a  GPS  time  stamp.  The  slaves  search  their  buffers  for  the  400  milli  second 
sample  that  corresponds  to  the  master's  time  stamped  FFT.  The  slaves  then  FFT  their 
own  data  and  correlate  these  samples  using  a  Complex  Ambiguity  Function  (CAP)  to 
determine  the  TDOA.  The  TDOA,  EDO  A,  GPS  position  of  the  observer,  velocity  of  the 
observer,  and  timestamp  are  returned  to  the  master.  The  master  collects  the  information 
and  sends  the  package  to  a  standard  Navy  TAC3/4  workstation  which  calculates  the 
geolocation.  A  detailed  reference  of  the  preliminary  hardware  configuration  and 
specifications  of  the  ARL:UT  system  can  be  found  in  Reference  20. 

Like  any  system  that  depends  on  the  Global  Positioning  System,  this  system  is  very 
vulnerable  if  GPS  were  lost  or  significantly  degraded.  While  relatively  reliable  sensor 
information  such  as  position  and  velocity  might  be  available  from  inertial  systems  on  the 
sensor's  platform,  the  need  for  a  reliable,  common  time  standard  between  sensors  so  that 
data  can  be  buffered  and  retrieved  at  common  times  would  still  be  required.  While  the 
vulnerabilities  of  GPS  is  an  issue  well  beyond  the  scope  of  this  discussion,  it  must  be 
considered  when  weighing  the  value,  limitations,  and  capabilities  of  such  a  system. 


m.  TDOA  AND  FDOA  MEASUREMENTS 

A.  BACKGROUND  AND  PURPOSE 

The  main  purpose  of  this  chapter  is  to  describe  how  the  TDOA  and  FDOA 
measurements  are  made  and  the  assumptions  behind  those  measurements.  Much  of  the 
next  section  is  a  basic  summary  of  several  signal  processing  theories  and  methods.  If  the 
reader  is  not  interested  in  these  details,  the  first  three  sections  of  this  chapter  may  be 
omitted  with  the  following  in  mind.  The  main  importance  of  this  chapter  is  to  present  the 
measurement  of  TDOA  's  and  FDOA  's  from  a  stochastic  RF  reception  as  a  statistical 
process  and  to  get  an  understanding  of  the  assumptions  required  to  formulate  these 
processes.  The  measurement  errors  in  any  TDOA  or  FDOA  are  multivariate  functions  of 
many  parameters — many  of  which  are  random  variables  themselves.  Measurements  are 
further  confounded  by  the  sensitivity  and  accuracy  of  the  method  employed  to  determine 
the  TDOA  or  FDOA  measurement. 

JBSTAFSim  does  not  assume  any  one  particular  method  to  measure  the  TDOA  or 
FDOA;  nor  does  JBSTAFSim  assume  any  properties  with  respect  to  the  signal,  the  noise 
in  the  signal,  or  the  noise  in  the  environment — the  random  variables  which  have  been 
shown  in  the  works  listed  below  to  function  as  parameters  which  describe  the 
measurement  error  estimate.  Rather,  JBSTAFSim  allows  the  user  to  set  the  measurement 
error  variance  estimate  and  correlation  of  the  measurements.  Since  these  measurements 
are  modeled  statistically,  it  is  reasonable  to  assume  that  regardless  of  the  method  used  to 
obtain  them  or  the  many  parameters  which  might  describe  them,  the  error  in  their 
measurement  can  also  be  described  statistically.  By  doing  this,  JBSTAFSim  can  simulate 
these  measurements  and  indicate  their  effect  on  the  geolocation  process  with  minimal 
information  provided  with  respect  to:  the  system  performing  the  measurements,  the  signal 
of  interest,  and  the  many  random  variables  describing  the  received  signal  and  its 
environment.  Viewed  this  way,  the  model  maintains  integrity  while  also  allowing  it  to 
serve  in  the  most  general  sense. 

The  following  discussion  is  based  on  Reference  1,  an  NPS  thesis  by  LT  David  A. 
Streight  who  applies  Cyclostationary  techniques  to  measure  the  TDOA  of  an  emitter.  In 
describing  the  advantages  of  this  particular  method  LT  Streight  provides  an  excellent 
overview  of  the  basic  statistical  properties  and  assumptions  of  the  signal  involved  in 


determining  the  TDOA  and   FDOA  measurement  using  traditional   methods.      More 
information  about  these  process  can  be  found  in  detail  in  References  4,  5,  and  6. 


B.  BASIC  MODEL  AND  THE  GENERALIZED  CROSS  CORRELATION 

Most  TDOA  models  assume  a  stationary  signal  which  is  received  by  two  sensors  in 
the  presence  of  white  Gaussian  noise.  These  signals  have  been  mathematically  modeled  as 
a  function  in  the  time  domain: 

x(t)  =  s(t)  +  n1(t)  2-1 

y(t)=As(t-D)  +  n2(t)  2-2 

where  s(t)  is  assumed  to  be  statistically  independent  and  uncorrelated  with  independent 
noises  «,(/)  and  n2(t) .    A  is  the  complex  valued  relative  magnitude  and  phase  mismatch 

between  the  receivers,    D    represents  the  TDOA  between  the  signals.  [Ref    1]     The 
autocorrelation  and  cross-correlation  functions  are  given  by: 

^(r)  =  ^(r)  +  /2;(r)  2-3 

^,(r)  =  M2/tfr)  +  /^(r)  2-4 

Ryx(T)=ARs(T-D)  +  RnjT)  2-5 

and  the  spectral  density  functions  are: 

Sx(f)  =  Ss(f)  +  Sni{f)  2-6 

S,(/H^(/)  +  ^(/)  2-7 

Sjf)  =  A-Ss(f)-e-J2*D+SniJf)  2-8 
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Using  these  relationships  and  assumptions,  the  parameter  D,  the  TDOA,  can  be 

obtained  by  using  correlation  techniques.    The  most  basic  technique  is  the  Generalized 
Cross  Correlation  method.    First,  note  that  (by  2-5)  R^ir)  peaks  when  t-D .    Next 

assuming  that  /  of  the  received  signal  is  within  a  finite  Bandwidth  B  around  the  carrier 

frequency  f0,  the  ratio  of  S^f)  to  Sx(f)  is: 

Sjf±_\A.e-™      /o-4</^/0+#  2-9 


sSf) 


2     "      '"     2 
otherwise 


The  inverse  Fourier  transform  of  this  ratio  is  equal  to: 


/o-S/2   WJ 

which  peaks  at  r  =  D  and  provides  the  TDOA  measurement.  [Ref.  1] 

C.  COMPLEX  AMBIGUITY  FUNCTION 

The  Generalized  Cross  Correlation  method  as  described  above  works  well  in  the 
case  of  static  transmitters  and  receivers.  However,  as  shown  in  Reference  5,  in  order  to 
determine  the  TDOA  when  there  exists  a  Doppler  shift,  the  spectrum  of  one  of  the  signals 
must  be  translated  in  frequency  equal  to  the  FDOA  measured  between  the  observers.  This 
can  be  shown  by  rewriting  the  generalized  model  above  so  that: 

x(t)  =  s(t)  +  n,(t)  2-11 

y(t)=As(t-D)e-j2^t  +n2(t)  2-12 

where  fd  is  the  Doppler  difference  measured  between  the  two  observers.    This  requires 
the  maximization  over  the  two  parameters  D  and  fd  which  leads  to: 
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J 

A(DJd)  =  jx(t)/(t-D)e-Wdt 


2-13 


Notice  that  at  fd  -  0,  the  CAF  reduces  to  the  Generalized  Cross  Correlation  method.  In 
Reference  5,  Stein  shows,  given  these  assumptions,  the  variance  of  the  estimates  for  the 
TDOA  and  FDOA  can  be  related  to  the  signal  integration  time,  the  signal  bandwidth,  and 
the  noise  bandwidth: 


TDOA 


J_  1 


2-14 


a 


FDOA 


\ 1_ 


2-15 


where: 


B=  the  noise  bandwidth  at  the  input  of  both  receivers, 


A 


-co 

\Ws(f)df 


2-16 


with  Ws(f)  =  the  spectral  density  of  the  signal  as  shaped  by  the  receiver, 


Te  =rms.  integration  time  and, 


r"2 


i     i       i 

—  +  —  +  — 

7\      Yi      YxYi 
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1 
with  —  =  the  signal  to  noise  ratio  (SNR)  for  each  receiver  [Ref  1] 

Y, 
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With  a  little  sensitivity  analysis  of  the  above  equations,  it  can  be  seen  that  the  driving 
component  in  the  standard  deviation  of  the  TDOA  measurement  is  SNR.  In  an  idealized 
case,  below  lOdB  SNR  and  with  a  400  ms  integration  time  of  the  ARL:UT  system,  the 
standard  deviation  of  a  single  TDOA  measurement  exceeds  100  meters.  "With  three 
measurements  to  provide  a  fix,  all  within  this  idealized  case,  the  geolocation  would  have  a 
standard  deviation  of  173  meters"  [Ref  1],  a  large  best-case  error. 


D.  CYCLOSTATIONARY  DETERMINATION 

LT  Streight's  thesis  applies  Cyclostationary  techniques  to  measure  the  TDOA. 
This  is  a  robust  approach  which  is  able  to  perform  well  despite  both  low  SNR  and  the 
presence  signals  which  may  or  may  not  intentionally  seek  to  jam  the  observers'  receivers. 
By  taking  advantage  of  the  cyclic — time  periodic — nature  of  most  signals,  a  three 
dimensional  correlation  in  the  frequency  and  cyclic  frequency  domains  can  be  performed. 
Thus,  where  the  signal  of  interest  and  noise  may  have  appeared  overlapping  in  the 
traditional  models,  they  can  still  be  separated  by  a  cyclic  autocorrelation  function  which 
can  be  thought  of  as  a  special  autocorrelation  function  which  produces  spectral  lines  at 
frequencies  indicative  of  the  signal's  periodic  nature.  In  the  frequency  domain,  the  cyclic 
autocorrelation  function  is  replaced  by  the  spectral  correlation  function.  The  Spectral 
Correlation  Function  (SCF)  is  a  correlation  in  the  frequency  domain  which  peaks  at  the 
nominal  spectral  frequency  and  cyclic  frequency  in  what  is  called  a  bi-frequency  plane. 
See  Figure  2-1. 
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Figure  3-1  Signals  plotted  in  the  Bi-frequency  Plane  (From  Ref.  1) 

This  method  allows  for  dramatically  better  results  in  the  presence  of  greater 
interference  and,  as  seen  above,  can  provide  a  classification  of  the  modulation  type.  The 
model  of  the  signal  at  any  pair  of  receivers  remains  as  in  2-1  and  2-2, 

*(/)  =  5(/) +  «,(/) 

y{t)=As(t-D)  +  n2{t), 

and  assumes  the  same  statistical  assumptions  with  respect  to  the  properties  of  the  signal 
and  noise  except  that  now  the  noise  may  now  contain  co-channel  interference  which  could 
spectrally  mask  the  signal.  As  before,  auto  and  cross-correlation  functions  are  used  in 
conjunction  with  the  spectral  density  functions  to  measure  the  TDOA;  however,  now 
these  correlated  values  are  aligned  by  the  cyclic  frequency  characteristic  of  the  desired 
signal.  While  highly  accurate  in  the  presence  of  interference,  this  method  is 
computationally  expensive.  LT  Streight  employs  a  computationally  efficient 
cyclostaionary  TDOA  algorithm,  SPECCOA,  "which  may  not  be  as  be  as  accurate  as 
others  but  is  the  best  choice  for  the  efficiency  required  in  a  tactical  system. "[Ref  1]  While 
not  implemented  by  the  ARL:UT  test  system,  this  method  could  be  implemented  as  the 
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algorithm  of  choice  in  the  next  system  or  in  future  systems.  It  is  presented  here  in  detail, 
along  with  the  GCC  and  CAF,  to  point  out  the  multitude  of  processes,  parameters,  and 
assumptions  required  to  estimate  the  TDOA  or  FDOA  measurement  error  variance. 
Regardless  of  the  method  used  to  obtain  them,  the  TDOA  and  FDOA  measurement  errors 
can  be  represented  as  a  stochastic  process. 


E.  IMPLICATIONS  OF  SIGNAL  PROCESSING  TO  GEOLOCATION 

Throughout  this  chapter,  discussion  has  focused  on  the  signal  models  and  their 
assumptions.  As  noted  previously  Stein,  in  Reference  5,  gives  the  TDOA  and  FDOA 
estimate  measurement  errors  based  on  a  few  parameters  and  several  assumed  conditions. 
Dr.  Petre  Rusu  and  Dr.  Lisa  Giulianelli  of  ARL:UT,  in  References  7  and  8,  follow  the 
work  of  Azaria  and  Hertz,  in  Reference  4,  to  find  the  solution  for  the  correlation  between 
TDOA  measurements  as  a  function  of  SNR,  observation  interval  length,  and  noise  and 
signal  bandwidths.  Unfortunately,  many  of  these  parameters  that  go  into  these 
formulations  are  in  fact  stochastic  themselves.  Further,  in  Reference  7,  the  enabling 
assumptions  with  regard  to  the  estimation  of  the  measurement  error  correlation  begin  with 
those  of  equations  2-1  and  2-2,  namely,  that  the  signal  received  at  any  observer  can  be 
described  as  the  summation  of  two  parts — a  deterministic  signal  plus  Gaussian  noise.  The 
noises  are  considered  ''uncorrelated  with  the  signal  and  with  each  other  across 
observers"[Ref  7].  While  these  assumptions  like  these  are  critical  for  the  necessary 
mathematical  formulations  like  those  in  sections  B  through  D  above  and  in  works  such  as 
Reference  7,  in  an  operating  environment,  these  models  will  fail  to  represent  the  one 
stochastic,  and  most  important,  variable  present  at  the  observers'  receivers — the  actual, 
"real  world"  signal  itself.  Viewing  each  signal  received  at  each  of  the  observers  as  a 
stochastic  process,  the  TDOA  and  FDOA  measurement  process  should  not  be  thought  of 
as  that  which  can  parameterize  the  measurement  errors;  rather,  it  is  the  randomness 
nascent  in  each  observer's  received  signal  that  introduces  the  measurement  errors. 

Since  quantifying  the  impacts  of  the  estimate  covariance  and  correlation  of  TDOA 
and  FDOA  measurement  errors  on  the  geolocation  process  is  the  goal  of  these  analyses, 
this  discussion  might  seem  to  be  disputing  a  fine  point.  However,  it  demonstrates  the 
evaluative  value  of  JBSTAFSim  as  a  tool  for  assessing  the  impact  of  any  assumption  or 
model  of  the  TDOA  and  FDOA  measurement  errors  on  the  geolocation  process. 
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JBSTAFSim,  then,  functions  as  an  "If — Then"  tool.  The  "If — Then"  portion  of  the 
simulation  is  remarkably  simple:  given  a  network  and  allocation  of  observers,  if 
measurement  errors  are  introduced  with  a  given  probabilistic  description,  then,  the 
expected  Geolocation  results  are  as  follows.  JBSTAFSim,  then,  exploits  one  of  the 
greatest  advantages  of  simulation,  rather  than  redefine  our  mathematical  models  or 
simulate  only  the  equations  of  our  models,  let  us  instead  violate  our  own  assumptions, 
stochastically,  and  examine  the  results  within  an  empirical  simulation  space.  The  results  of 
such  an  investigation  are  useful  for  anyone  who:  is  designing  a  system  such  as  this  and 
wishes  to  evaluate  or  improve  its  effectiveness,  using  a  system  such  as  this  for  tactical 
decisions  and  wishes  to  know  its  capabilities  and  limitations,  and  for  anyone  sponsoring 
the  development  of  such  a  system  and  wishes  to  know  how  added  capabilities  or 
improvements  might  improve  the  warrior's  use  of  such  a  system. 
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IV.  EMITTER  GEOLOCATION 

A.  BACKGROUND 

The  TDOA  and  FDOA  measurements  computed  by  the  signal  processing  methods 
described  in  the  previous  chapter  form  the  basis  for  the  measurement  data  required  to 
solve  for  the  emitter's  location.  These  quantities  measured  in  seconds  and  in  Hertz  can  be 
represented  formally  in  terms  of  length  (meters)  and  velocity  (meters/second)  by  the 
following  algebraic  equations: 

TDOAff=<{tt-tJ) 

tdoa, = v(*, -*.)* +u -y.r +(*. - r -fa-*.)2 +h -*r +b -j  4_1 


FDOA,  =  -/, 


>: 


Vj 


Where  the  subscripts  /  and  j  represent  the  observer  pair  forming  the  differences  and  s 
denotes  the  signal's  source.  An  observer's  three  dimensional  location  in  space  is  defined 
by  (x,y,z),  and  r  and  v  represent  the  radial  and  velocity  vectors  of  the  target  or 
observers.  Therefore,  the  TDOA  and  FDOA  can  be  thought  of  as  spatial  differences  and 
relative  velocity  differences,  rather  than  in  time  and  frequency  units 

Much  has  been  written  and  theorized  about  the  best  way  to  solve  for  the  state  of 
an  emitter  given  TDOA  information.  While  the  TDOA  measurement  is  non-linear  in  the 
emitter's  state,  solutions  can  be  formulated  to  directly  solve  for  the  location  of  the  target 
as  linear  models  which  satisfy  non-linear  constraints.  These  solutions  usually  involve 
squaring  the  TOA  equation  which  results  in  two  locations  which  satisfy  the  TDOA  data. 
References  9,  10,  and  12  all  show  different,  but  mathematically  equivalent  algebraic 
methods  to  directly  solve  for  the  state  of  an  emitter.  In  particular,  in  Reference  10,  Dr. 
Petre  Rusu,  of  ARL:UT,  goes  through  great  work  to  linearly  propagate  the  estimate 
errors  through  his  solution  to  provide  an  uncertainty  estimate  for  his  solution  to  the 
emitter's  state.  In  Reference  11,  Dr.  Rusu  further  proves,  given  consistent  weighting 
matrices,  the  theoretical   equivalence  of  the  TOA  direct   solution   and   least   squares 
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formulation.  References  9  and  13  are  specifically  geared  to  the  algebraic  solution  of  GPS 
equations;  however,  their  techniques  can  be  applied  to  any  TDOA  geolocation  problem. 

In  JBSTAFSim,  the  method  of  Reference  13,  by  Chaffee  and  Abel,  is  used  to 
provide  an  initial  starting  position  for  the  linearized  least  squares  approach.  This  method 
is  very  elegant  and  numerically  robust.  Further,  it  provides  a  useful  understanding  of  the 
solution  ambiguity  by  not  immediately  squaring  the  TOA  equation.  Thus,  this  method  was 
chosen  for  its  numerical  superiority  and  for  its  elegant  resolution  of  the  ambiguity  of 
emitter  location  The  linearized  least  squares  solution  uses  this  direct  solution  as  its 
starting  point  from  which  it  seeks  a  solution  to  minimize  the  sum  of  squared  errors.  After 
a  fix  is  computed,  the  solver  then  uses  the  previous  solution  as  the  starting  position  for  the 
subsequent  solutions  until  a  request  is  made  for  another  direct  solution  to  update  the 
current  solution.  The  periodicity  of  these  updates  may  be  set  by  the  user.  In  an  actual 
geo-location  system,  direct  solution  updates  might  be  requested  as  a  function  of  many 
complex  variables  and  parameters,  for  instance:  quality  of  the  previous  fix,  the 
dimensionless  shock,  time  since  the  last  fix,  by  request  of  the  system  user,  or  quality  of  the 
received  signal.  JBSTAFSim  seeks  to  minimize  the  sum  of  squared  errors  by 
implementing  a  Square  Root  Information  Filter  (SRIF).  The  SRIF  was  chosen  for  its 
numerical  superiority  and  its  ability  easily  include  a  priori  statistics  and  function  as  a 
numerically  stable  extended  Kalman  filter.  The  SRIF  in  JBSTAFSim  can  perform  either  as 
a  single  epoch  state  solver  or  as  an  extended  Kalman  Filter,  "batch",  processor. 
JBSTAFSim's  SRIF  does  not,  however,  incorporate  a  movement  model  that  considers  a 
target  which  moves  with  time.  While  a  movement  model  would  no  doubt  provide  the 
most  robust  approach,  given  limited  communications  between  observers,  the  complexity 
of  such  a  model,  the  non-linearity  of  the  measurements,  and  the  purpose  of  the  system  to 
act  as  a  locator  and  not  as  a  fire  control  system,  the  batch  process  mode  without  a 
movement  model  should  be  considered  adequate  for  those  targets  with  zero,  or  near  zero, 
velocities.  The  development  of  movement  models  for  particular  classes  of  targets  would 
be  a  large  undertaking,  however,  their  design  could  be  facilitated  through  empirical 
simulation  by  extending  the  JBSTAFSim  solver. 

Finally,  and  as  will  be  discussed  in  the  next  chapter,  JBSTAFSim's  simulation 
results  should  not  be  viewed  as  a  measure  of  the  ARLUT  system  In  particular, 
ARLUT's  system  is  currently  in  the  test  and  development  phase  in  which  its  Engineers 
and  Scientists  are  considering  a  multitude  of  design  issues  like  the  ones  presented  in  this 
thesis.     Therefore,  JBSTAFSim's  results  are  not  necessarily  identical  with  that  of  the 
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the  exact  solution  for  the  GPS  pseudorange  equations.    This  method  provides  a  greater 
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unuerSianuing  Or  me  geometry  inVOiVeu  m  trte  lvkj/\  equations  anu  can  oe  applied  in  a 
numerically  stable  fashion.  Dr.  Tolman  of  ARL:UT  provides  a  simple  explanation  of 
Chaffee  and  Abel's  formulation  in  Reference  14.  This  section  provides  a  quick  summary 
of  this  direct  solution  method  based  on  both  of  these  references. 

Following  the  notation  and  example  of  Dr.  Tolman's  notes  in  Reference  14,  let  the 
emitter  be  located  at  position  f  and  define  an  observer  to  be  located  at  position  ft .   The 

range  between  one  observer  and  the  emitter  is  defined  by: 

/Hfc-i  4-3 

The  TOA  equation  is  defined  as: 

7O4=c/I.=||^-r|  +  cr0=^.+rf0  4-4 

where  *0  is  the  unknown  time  of  transmission  of  the  signal.  This  equation  holds  for  any 
change  in  the  origin  of  time.  The  TDOA  is  simply  the  difference  of  two  TOA's: 

Algorithms  in  References  9  and  1 0  begin  by  squaring  the  TOA  equation  above.  Chaffee 
and  Abei  do  not.  Instead,  iabei  one  receiver  I,  the  master.  Let  aii  other  receivers  be 
denoted  with  /  s  1, ...  M — 1 . 

ctt  -ctj  -(ct0 -<*,)=  |/J  -fj  -(r-rj|  4-6 
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Further,  define: 


A  =  ct  -ct, 


ri  =ri~ri 


r  -r  -r, 


This  makes  the  TOA  equation: 


A.+lhHK'-'ll 


4-7 


r.ll2 


Square  this  equation  and  note  \\r'\\    cancels. 


ir1|2+2rV=A'+2|HA, 


4-8 


Stack  this  equation  for  M-l  differences  yields: 


f'    A. 


%'•%'- A* 


4-9 


For  continuity,  Chaffee  and  Abel  describe  this  set  of  equations  as: 


A. 


u 

lli/ll 


3 

A.. 


S, 

A. 


4-10 


Where  (v,w)  represents  the  Lorentz  inner  product  on  9?4  defined  as:    if  v  =  (x' ,a)  and 


w 


=  [y',c)   are  vectors,  then  (v,w)  -  x'y-ac.     In  addition,  also  note  Z 


S,  =f',and  u  =  F.\Ref.  13] 

Note  that  these  two  equations,  4-9  and  4-10,  are  both  of  the  form  AX  =  B  .  Thus, 
this  appears  to  be  a  linear  problem  in  X,  but  the  emitter  location  solved  by  TDOA 
techniques  is  not  linear.   In  this  case,  the  fourth  component  of  X  is  non-linear;  this  non- 
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linearity  will  act  as  a  separate  constraint  on  the  solution  of  the  linear  problem.  In  the 
overdetermined  case  when  there  are  four  equations  (M=5  observers),  the  non-linear 
constraint  is  automatically  satisfied,  however,  we  rarely  expect  to  be  in  an  overdetermined 
case.  Therefore,  the  matrix  A  is  singular,  but  still  can  be  inverted  as  a  generalized 
inverse:     A~z  -  At(AAt)    .    The  best  way  to  obtain  A~g ,  however,  is  by  the  singular 

value  decomposition  (SVD)  algorithm.  More  compelling,  the  SVD  algorithm  directly 
provides  the  null  space,  a1 ,  required  to  satisfy  the  non-linear  constraint.  Since  A  is  3x4, 
and  we  assume  rank  3,  then  the  null  space  of  A  is  one  dimensional  and  the  complete 
solution  for  X  is: 

X  =  A~gB  +  AaL  4-11 

Where  aL  is  a  unit  vector  spanning  the  null  space  of  A  and  1  is  a  real  variable.  To 
prove  this  is  the  solution  to  AX  =  B ,  multiply  both  sides  of  the  solution  by  A  and  recall 
Aa1  =0, 

AX  =  AA1 \AATY B  +  XAaL  =  AX  =  B  4-12 

Given  the  above  equation,  it  is  evident  that  the  constant  /  is  arbitrary  with  respect 
to  the  linear  equation;  however,  from  equation  4-9,  the  non-linear  constraint  requires  the 
fourth  component  of  X  be  equal  to  the  magnitude  of  the  vector  of  its  first  three 
components:   ||r'|  -(X4)   =0. 

Now,  consider  the  first  three  dimensions  or  components  of  the  solution  for  X  — 
Chaffee  and  Abel  call  this  vector:  uA  -  A'yAA')    u  -  A~gB .    The  solution  to  the  linear 

equation  has  given  us  this  vector  which  describes  the  emitter's  location  uniquely  in  a 
three-dimensional  subspace  orthogonal  to  a  ,  but  cannot  provide  any  information  to  the 
ambiguous  component  normal  to  this  subspace,  i.e.  along  a~ .  The  non-linear  constraint 
must  be  used  to  determine  this  component  defined  by  X .  Substituting  the  four 
dimensional  solution  to  the  linear  equation,  Chaffee  and  Abel  label  wA  ,  plus  aci1  back 

into  the  Lorentz  functional,  yields  a  quadratic  equation  for/. : 

(a1  ,al)X2  +2Z(wA,a1)  +  (wA,wA)  =  0  4-13 
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which  provides  a  direct  solution  for  A  . 

Thus,  to  summarize  the  result  and  combining  the  two  notations: 


=  wA  +  Aa    +  Z 


4-14 


wA  +  Aa 


The  quantity  b  in  Equation  4-14  is  called  the  bias  which,  used  in  GPS 


where     [wA  +  Aa1  J  >  0     which    is    from    the    constraint    on    the    solution    vector 
P~ 

IMI 

equations,  represents  the  offset  of  time  between  the  master's  and  emitter's  time  reference. 
The  vector  f  gives  the  emitter's  location. 

Like  all  direct  solutions,  the  impact  of  squaring  the  TOA  equation  results  in  two 
possible  solutions  for  the  emitter  location.  For  GPS,  this  is  not  a  problem  since  one 
solution  is  on  the  face  of  the  earth  and  the  other  lies  opposite  to  the  satellite  plane  in 
space.  In  a  tactical  environment,  however,  this  problem  must  be  solved  with  some  prior 
knowledge  of  the  target  of  interest  or  with  some  common  sense  with  respect  to  the 
components  which  describe  the  target's  location. 

C.  LINEARIZED  LEAST  SQUARES  APPROACH 

The  linearized  least  squares  approach  to  the  geolocation  of  an  emitter  provides  a 
statistically  robust  estimate  of  the  emitter's  state  and  covariance  uncertainty  estimate. 
While  the  method  does  require  an  initial  guess  with  respect  to  the  location  of  the  emitter, 
in  many  cases  this  can  be  provided  by  picking  a  point  in  the  area  of  observers  or  by  first 
solving  for  the  state  via  one  of  the  direct  methods  mentioned  above.  Certainly,  the 
performance  of  linear  least  squares  algorithms  is  well-known  and  stable  given  that  care  is 
taken  into  their  formulation.  In  this  regard,  JBSTAFSim  utilizes  a  Square  Root 
Information  Filter  (SRIF),  see  Appendix  A  for  a  detailed  description,  which  has  been 
proven  to  be  more  numerically  stable  and  accurate  than  other  methods.  In  addition,  the 
SRIF  allows  for  elegantly  simple  inclusion  of  prior  estimates  to  be  included  with  the 
current  state  estimate  for  batch  processing  as  an  extended  Kalman  Filter.  Further  a  hybrid 
linearized  least  squares  solution  aided  by  a  direct  solution  will  most  likely  be  the 
geolocation  method  of  choice  by  the  ARL:UT  system.    In  the  mixed  TDOA/FDOA  case, 
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JBSTAFSim  also  enhances  the  numerical  accuracy  of  the  SRIF  solution  by  scaling  the 
partials  matrix  and  TDOA  and  FDOA  data  so  that  every  equation  is  dimensionless  and  on 
unity  order.  Appendix  B  contains  a  detailed  description  of  the  scaling  of  mixed 
TDOA/FDOA  data  in  order  to  improve  the  numerical  accuracy  and  stability  of  a  linearized 
least  squares  solution. 

Consider  in  detail  the  mixed  TDOA/FDOA  problem  as  described  in  Reference  19: 

Observation  equations: 


c-mOAlJ(F)  =  R(')-R 
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Measurement  covariance  matrix: 
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Linearized  least  squares  equation: 
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Direction  cosines  are  denoted  by  a ,  and  R  is  the  observer-target  separation.  In 
the  state  vector,  dx  and  dv  refer  to  the  offset  from  the  position  and  velocity,  respectively, 
of  the  nominal  state,  and  the  subscripts  p   and  _L   refer  to  components  parallel  and 

perpendicular,  respectively,  to  the  line  of  sight.  Thus,  for  mixed  T/FDOA,  the  above 
partials  matrix  is  2(M  -  l)xN ,  where  M  equals  the  number  of  observers  and  N  is  the  size 
of  the  unknown  emitter's  state.  When  only  TDOA  information  is  used  for  geolocation, 
the  partials  matrix  is  reduced  to  (M  -  l)xN  ,  where  N  <  3 ,  and,  thus,  no  information  with 
respect  to  the  target's  velocity  can  be  given.  For  the  purposes  of  this  thesis,  only  cases 
which  are  evenly  determined  are  considered.     The  overdetermined  solutions  will  most 
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likely  always  dominate  the  evenly  determined  ones;  further,  this  method  of  analysis 
assumes  that  resources  are  typically  scarce  and,  thus,  the  overdetermined  case  event  is 
unlikely  to  occur.  With  this  philosophy  in  mind,  JBSTAFSim  is  designed  to  produce  two- 
or  three-dimensional  fixes  in  the  TDOA  only  case  and  four-  or  six-dimensional  fixes  in  the 
mixed  T/FDOA  mode. 

In  addition  to  the  solution  to  the  state  of  the  target,  the  linearized  least  squares 
approach  also  provides  an  estimate  of  the  error  co variance  of  the  state  estimate.  Viewing 
4-18  as  the  set  of  equations  Ax  =  b ,  the  estimate  of  the  state  vector  takes  on  the  form 
{ATM~]A)    -P.    If  we  are  interested  in  only  the  first  n  components  of  the  solution 

vector  and  take  the  first  men  components  of  P  correspondingly;  then,  given  the  implicit 
assumptions  of  Normality  in  the  least  squares  and  extended  Kalman  filter,  it  can  be  shown 
that: 

{*.-mTpA**-m)~£  4"19 

where   ju    represents  the  true  location  of  the  target.    Since  this  co  variance  matrix 

parameterized  by  n  is  always  positive  definite,  it  can  be  said  that  its  column  vectors  form 
the  basis  of  the  eigenspace  Pn  and  the  corresponding  eigenvalues  are:    av..an.    Thus  a 

hyperellipsoid  defined  by: 
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forms  surfaces  of  equal  probability  densities.  Therefore,  for  a  certain  £ ,  the  probability 
that  a  point  will  he  inside  the  ellipsoid  can  be  computed.  For  example,  if  an  n=2 
dimensional  fix  with  an  £  =  3  sigma  error  ellipse  is  desired,  the  ellipse  defined  by  the 
covariance  matrix  would  yield  an  ellipse  for  which  there  would  be  a  99.7%  probability  that 
the  target  lies  inside  the  ellipse.  Thus,  the  size  of  the  ellipsoid  created  by  the  covariance 
matrix  can  be  thought  of  as  a  measure  of  the  uncertainty  of  the  solution. 

In  addition  to  generating  error  ellipsoids,  the  covariance  matrix  also  provides  a 
means  to  measure  the  amount  by  which  the  observer  geometry  is  degrading  the  fix.  This 
measurement  is  often  referred  to  as  the  geometric  dilution  of  precision  (GDOP).  Again, 
given  that  the  residuals  of  the  state  vector  are  Normally  distributed,  the  covariance  matrix 
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is  a  measure  of  their  statistical  properties.  If  these  errors  are  small,  a  linear  approximation 
which  relates  them  to  the  observations  may  be  considered  satisfactory.  "Moreover,  if  the 
observation  residuals  are  Normally  distributed,  the  fix  residual  residuals  computed  via  the 
linear  approximation  are  automatically  Normally  distributed."[Ref.  15]  In  this  case  the 
linear  error  propagation  may  be  defined  from  Reference  1 5  as: 

Cy  =  AMAT.  4_2i 

The  total  error  on  the  fix  as  the  root-mean-square  deviation  is 


^Tr(Cy)  =  ^Tr(AMAT)  4-22 

which  provides  the  "total  fix  position  error."  Dividing  this  quantity,  as  suggested  above, 
by  the  total  measurement  error  provides  an  estimate  of  the  relative  dilution  of  precision: 

^Tr(AMAT)  4-23 

■y/MM) 

While  there  is  nothing  geometric  about  4-23,  it  places  all  the  geometric  errors  in  the 
\AAT)   factor  separately  from  the  stochastic  contribution  of  C   .     The  relationship 

between  GDOP,  measurement  error,  and  fix  error  has  been  described  as: 

G  position  *  GDOP  ■  CTmeaSurement  4-24 

Therefore,  the  GDOP  is  not  truly  geometric,  but  depends  on  the  stochastic  portion  of  the 
model.  It  is,  however,  "accurate  for  quantitative  predictions  within  a  class  of  stochastic 
models  defined  by  covariance  matrices  that  are  not  too  different  from  each 
other"[Ref.  15].  Further,  it  appears  that  based  on  the  numerical  experience  of  the 
ARL:UT  test  system,  that  it  is  a  good  tool  for  quantitative  prediction.  [Ref.  20]  As  will  be 
shown  in  Chapter  VII,  while  GDOP  provides  a  good  relative  measure  of  the  geometry,  a 
better  means  to  measure  the  degree  to  which  the  geometry  of  the  observers  is  needed. 
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V.  JAVA-BASED  STOCHASTIC  TDOA  AND  FDOA  SIMULATION 

(JBSTAFSEVf) 


A.  BACKGROUND 

JBSTAFSim  was  written  to  stochastically  explore  the  differences  and  similarities 
of  TDOA  and  TDOA/FDOA  GPS-assisted  geo-location  solutions.  It  assumes  a  Multi- 
platform,  GPS-assisted  architecture  similar  to  that  being  developed  by  ARL:UT  as 
described  in  the  previous  chapters.  JBSTAFSim  allows  the  user  to:  create  scenarios  of 
various  types  of  sensors  and  targets  defined  by  straight-forward  parameters;  distribute 
these  platforms  in  a  simulated  operating  environment;  dictate  solver  characteristics  for  the 
linear-least  squares  formulation;  and  adjust  the  measurement  error  distribution  and 
correlation.  While  the  purpose  of  JBSTAFSim  was  to  mainly  explore  the  differences 
between  the  two  types  of  solutions  in  the  evenly  determined  case,  it  could  also  be  used  a 
tool  to  investigate  and  experiment  with  various  solver  properties,  solution  methods,  or 
algorithms  in  a  variety  of  cases.  Further,  it  could  simply  be  used  to  experiment  with 
various  platform  combinations,  characteristics,  or  allocations  to  provide  an  understanding 
of  the  system's  strengths,  weaknesses,  and  limitations.  Finally,  JBSTAFSim  is  an  object- 
oriented  program  which  has  been  written  to  be  extendible.  Given  an  understanding  of  the 
model  as  described  in  the  following  sections,  the  simulation  can  be  extended  to  satisfy  any 
of  the  areas  of  study  discussed  above  or  modified  to  include  more  complex  sensor  or 
target  models. 

B.  DISCRETE  EVENT  SIMULATION  IN  JBSTAFSIM 

If  geo-location  of  an  RF  transmitter  using  TDOA  and  FDOA  techniques  was  a 
straight-forward  analytical  problem,  then  a  simple  closed-form  solution  could  be  produced 
on  paper.  Yet,  the  signal  processing  techniques  and  assumptions  presented  in  Chapter  III, 
show  that  any  realistic  geo-location  system  is  a  "system"  built  on  highly  complex  models 
so  complex  that  any  one  analytic  model  designed  to  mathematically  represent  the  whole 
system  will  fail  to  realistically  describe  the  real-world  system.  Simulation,  on  the  other 
hand,  provides  a  means  to  "evaluate  a  model  numerically,  and  . .  [gather  data]  to  estimate 
the  desired  true  characteristics  of  the  moder  [Ref  23].  Therefore,  simulation  provides  a 
means  to  examine  the  output  measures  of  effectiveness  by  numerically  exercising  the 
model  wfth  inputs  that  would  violate  the  assumptions  of  the  analytic  model  and,  thus, 
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provide  estimates  that  are  more  characteristic  of  the  real-world  system.  Given  this 
complex  mathematical  model  which  acts  dynamically  and  is  influenced  by  stochastic 
components,  JBSTAFSim,  then,  provides  an  evaluative,  numeric  tool  which  can  evaluate 
and/or  estimate  the  performance  of  a  Multi-platform  GPS-assisted  geo-location  system  as 
it  would  perform  in  a  realistic  environment — an  analysis  well  outside  and  beyond  the 
bounds  of  that  provided  by  any  mathematical  formula,  theorem,  or  set  of  equations. 

In  Discrete  Event  Simulation  (DES),  simulated  time  only  passes  when  events  are 
not  occurring.  "In  mathematical  terms,  we  might  say  that  the  system  can  change  [state]  at 
only  a  countable  number  of  points  in  time"  [Ref  23].  Events  are  listed  on  a  master 
schedule  in  time  and  priority  order.  The  entity  which  passes  simulated  time  and  ensures 
that  events  are  executed  in  simulated  time  is  often  referred  to  as  the  Time  Master.  When 
no  events  are  occurring,  the  Time  Master  simply  advances  the  simulation  clock  to  the 
scheduled  time  for  the  next  event  on  the  event  list.  The  Time  Master  then  allows  that 
event  to  execute,  this  event  may  cause  other  events  to  execute  or  be  scheduled  with  the 
Time  Master.  When  the  event  is  complete,  the  Time  Master  takes  control  again  and 
executes  the  next  event  on  the  event  list,  advancing  time  as  appropriate.  When  there  are 
no  more  events  left  to  schedule  or  execute,  the  simulation  is  finished. 

While  Java  is  not  strictly  a  simulation  language,  as  Java  gains  popularity,  new 
packages  are  being  developed  for  many  Java-based  applications,  including  simulation. 
JBSTAFSim  utilizes  a  simulation  package  named  SIMKIT  which  was  developed  at  Naval 
Postgraduate  School  by  LT  Kirk  Stork,  USN,  in  conjunction  with  Professor  Arnold  H. 
Buss.  SIMKIT  is  a  package  designed  to  facilitate  Discrete  Event  Simulation  (DES).  The 
version  of  SIMKIT  used  in  JBSTAFSim  is  a  development  release  written  in  the  JDK  1.0.2 
and  utilizes  the  version  of  Object  Space's  Java  Generic  Library  (JGL)  written  for  the  JDK 
1.0.2.  Since  SIMKIT  is  under  development,  further  references  to  SIMKIT  shall  refer  to 
the  development  release  version  of  SIMKIT  used  in  JBSTAFSim.  Detailed  information 
about  the  JGL  and  SIMKIT  can  be  found  at  www.objectspace.com  and 
http://131. 120.  142. 115/~stork/simkit  home. html,  respectively. 

In  addition  to  providing  the  DES  mechanization,  SIMKIT  also  provides  a 
"RandomStream"  class  to  produce  psuedorandom  variables.  While  Java  does  implement  a 
psuedorandom  number  generator,  SIMKIT's  RandomStream  class  not  only  provides  more 
suitably  random  numbers,  but  also  provides  the  user  a  library  of  random  variate 
distributions.  JBSTAFSim  employs  several  of  these  psuedorandom  number  generators  for 
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platform  movement  and  makes  use  of  SEVQCIT's  random  uniform  [0,1]  generator  to 
produce  the  random  distributions  for  TDOA  and  FDOA  measurement  errors. 

C.  MODEL  ASSUMPTIONS 

Before  going  too  deeply  into  model  assumptions  and  sensor  models,  it  is  important 
to  understand  that  JBSTAFSim  was  designed  to  provide  a  stochastic  simulation  for  the 
comparison  of  TDOA  and  TDOA/FDOA  geolocation  systems.  Like  any  good  model,  the 
model  strives  for  simplicity  wherever  possible.  As  described  below,  platforms  which  move 
in  the  simulated  environment  are  fairly  simple  yet  retain  their  basic  qualities  which 
differentiate  them  from  each  other.  Further,  the  logic  of  the  solver  is  a  very  limited  to 
solving  for  a  target's  position  given  the  data  provided  by  the  sensors.  The  solver  expects 
an  evenly  determined  case  and  it  does  not  accept  user  guidance  during  the  simulation  nor 
does  it  make  decisions  using  a  complex  algorithm  or  neural  network.  Functions  such  as 
these  are  those  which  are  assumed  to  be  in  place  in  an  operational  system.  In  addition  to 
assuming  a  geo-location  system  which  is  properly  functioning  and  that  the  master  has 
good  communications  with  its  slave  sensors,  the  following  sections  describe  the  clarifying 
assumptions  made  in  JBSTAFSim. 

1.  GPS  Data 

First,  JBSTAFSim  assumes  that  the  platform  GPS  positions  and  velocities 
represent  truth.  For  the  model,  this  is  a  necessary  condition  since  the  model  must  establish 
some  representation  of  ground  truth.  In  reality,  this  is  also  a  fairly  good  assumption  given 
that  a  survey-grade  GPS  unit,  which  is  fully  Precise  Positioning  Service  (PPS)  and 
precision  (p-code)  capable,  is  implemented  in  the  ARL:UT  test  platform  With  this  type  of 
equipment,  GPS  fixes  are  considered  to  be  accurate  to  within  10  meters.  With  respect  to 
velocities,  the  manufacturer  of  the  GPS  receiver  in  the  ARL:UT  prototype  system  asserts 
velocity  measurements  accuracy  of  0.01  meters  per  second.  [Ref  24] 

As  a  consequence  of  using  GPS  information,  the  coordinate  system  used  by 
JBSTAFSim  and  the  ARL:UT  system  is  the  earth-centered-earth-fixed  coordinate  system 
as  established  by  the  World  Geodetic  System  of  1984  (WGS84  ellipsoid).  This  system 
places  all  platforms  in  an  ellipsoid  defined  by  three  coordinates  (X,  Y,  Z)  in  units  of 
meters.   In  order  to  assist  the  user  in  inputting  data  and  for  graphically  displaying  results, 
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JBSTAFSim  receives  and  displays  data  in  coordinates  Latitude,  Longitude,  and  Height. 
All  calculations  and  results,  however,  are  measured  in  meters  in  the  WGS84  ellipsoid. 

2.  Platform  Models 

The  platform  models  which  function  in  the  JBSTAFSim  simulated  environment 
require  some  clarifying  assumptions  for  each  platform's  movement  within  this  simulated 
space.  In  order  to  understand  how  the  JBSTAFSim  platform  models  function  in  the  DES 
environment,  Figures  5-1  through  5-4  depict  event  graphs  for  each  platform. 

a.         Ship  Model 

The  ship  model  retains  the  most  common  features  of  all  five  platform 
models.  Like  all  platform  models,  the  ship  keeps  track  of  its  current  position  and  velocity 
in  earth  centered  earth  fixed  coordinates.  The  event  graph  for  a  ship's  movement  is  shown 
in  Figure  5-1.  Upon  creation,  any  ship  in  JBSTAFSim  is  provided  by  the  user  an 
operating  box  in  which  to  steam.  Also  at  creation,  the  ship  chooses  independent  random 
uniform  X  and  Y  coordinates  in  this  box  as  its  destination.  The  ship  calculates  the  amount 
of  time  it  will  take  to  arrive  at  its  destination  based  on  its  current  speed  and  distance  to  the 
destination.  The  ship  then  schedules  an  event  to  create  a  new  destination  when  it  arrives 
at  the  current  destination.  In  addition  to  destination  events,  ships  schedule  their  own 
change  speed  events.  New  ship  speeds  are  selected  randomly  from  a  uniform  distribution 
defined  by  the  low  and  high  speed  bounds  provided  by  the  user.  Occurrences  of  speed 
change  events  are  modeled  as  events  with  an  exponential  distribution  with  input  parameter 
provided  by  the  user.  The  first  speed  change  event  is  scheduled  upon  initialization.  When 
a  speed  change  event  occurs,  in  addition  to  changing  the  speed  of  the  ship,  the  waiting 
change  destination  event  is  canceled  and  a  revised  arrival  time  is  computed  based  on  the 
new  speed  of  the  ship.  A  new  change  destination  event  is  then  scheduled  for  this  revised 
time.  Finally,  the  speed  change  event  schedules  a  new  speed  change  event  to  occur  at  a 
random  time  provided  by  the  exponential  distribution. 
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Figure  5-1.  Ship  Model  Event  Graph. 


b.  Ground  Unit  Model 

The  ground  unit  model  is  slightly  more  complex  than  the  ship  model.  The 
event  graph  for  the  ground  unit  model  is  provided  in  Figure  5-2.  Like  the  ship  model  the 
ground  unit  keeps  track  of  its  position  and  velocity  and  schedules  change  destination 
events  and  change  speed  events  The  difference  with  the  ground  unit  is  that  it  occasionally 
ceases  to  move,  or  does  what  is  called  a  "take  break  event."  Break  events  are  scheduled 
as  a  random  process  with  an  exponential  distribution  with  input  parameter  provided  by  the 
user.  The  length  of  such  take  break  events  are  also  modeled  randomly  by  an  exponential 
distribution  with  input  parameter  provided  by  the  user.  When  a  take  break  event  occurs, 
in  addition  to  zeroing  the  speed  of  the  ground  unit,  the  waiting  change  speed  event  and 
change  destination  event  are  canceled.  A  revised  change  destination  event  is  scheduled  for 
the  remaining  transit  time  to  the  current  destination  but  only  after  the  amount  of  break 
time  has  elapsed.  Similarly  a  change  speed  event  is  scheduled  in  an  amount  of  time 
provided  by  the  random  distribution  plus  the  amount  of  rest  time  for  the  break.  Finally, 
after  the  rest  period  has  elapsed,  a  take  break  event  schedules  another  take  break  event  in 
a  random  amount  of  time  provided  from  the  exponential  distribution. 
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Figure  5-2    Ground  Unit  Model  Event  Graph. 


c  Aircraft  Model 

The  aircraft  model  is  more  simple.  Each  aircraft  patrols  a  circular  "orbit" 
defined  by  the  circle's  center  and  radius.  In  addition  to  these  parameters,  the  user  also 
defines  the  aircraft's  base  height  and  how  much  the  aircraft  deviates  from  this  base.  Since 
the  aircraft  moves  around  its  center  in  a  fixed  two-dimensional  circle,  the  aircraft's  speed 
is  vital  to  determining  where  on  the  circle  the  aircraft  lies.  For  any  circle  with  radius  r  the 

angle  6  at  the  center  of  the  circle  swept  by  an  aircraft  traveling  along  an  arc  length  s  is 

s 
simply  defined  as  0  =  — .    The  rate  at  which  this  angle  changes  can  be  simply  defined  by 


dO     1  ds 
the  first  derivative:     — -  =  — - 

dt      r  dt 


ds 
co .     Since  —  is  a  measure  of  the  magnitude  of  the 
dt 


aircraft's  velocity,  |v|,  co  can  be  computed  as  — .    With  this  value  of  co ,  the  aircraft's 

position  on  the  circle  can  be  computed  as  a  function  of  time.  Thus,  the  aircraft  model 
need  only  keep  track  of  co  and  update  its  value  when  |v|  changes.  In  order  to  prevent  two 
aircraft  from  having  the  same  or  near  similar  values  for  6 ,  each  aircraft  is  given  its  own 
offset  angle  chosen  randomly  over  the  interval  [0,  211].   With  this  model,  an  aircraft  need 
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only  update  its  position  on  its  circular  track  when  polled  and  then  pick  a  random  height 
adjustment  to  make  to  its  base  height.  However,  when  polled,  the  aircraft  must  also 
update  its  current  velocity  state.  In  order  to  compute  the  current  velocity  of  an  aircraft, 
the  aircraft  need  only  look  a  very  short  amount  of  time  into  the  future  and  draw  a  straight 
line  to  this  point  to  set  its  velocity  vector.  With  this  construct  in  mind,  the  event  graph 
below  is  simple  to  understand.  The  aircraft  need  only  update  positions  and  velocity  when 
polled  or  when  a  speed  change  event  changes  |vj .  As  before,  speed  change  events  occur 
randomly;  are  modeled  with  an  exponential  distribution;  and  cause  a  new  speed  change 
event  to  be  scheduled. 


Poll 
Aircraft 

/        Aircraft 
Creation 

Update       \ 

Position          J 

/v\. 

v     7 

Speed       \ 
Change         j 
Evert         / 

7 

Figure  5-3.  Aircraft  Model  Event  Graph. 


d  Satellite  Model 


While  much  could  be  done  to  model  satellites  of  different  orbital 
characteristics,  the  satellite  model  in  JBSTAFSim  is  remarkably  simple.  The  reason  for 
this  simplicity  is  for  the  model  to  serve  its  purpose  as  an  evaluative  simulation.  Certainly, 
more  complex  orbits  could  be  included  by  extending  the  current  model.  JBSTAFSim's 
satellite  model  assumes  a  satellite,  or  satellites,  in  a  circular  orbit  with  no  gap  in  coverage. 
The  model  only  requires  the  orbit  height,  starting  and  ending  footprints  of  the  region  of 
coverage.  Simulating  continuous  coverage,  the  satellite  flies  starting  point  to  endpoint  and 
returns  back  to  the  starting  point.    Velocity  information  need  not  be  entered  since,  for 


satellites  in  a  circular  orbit,  the  magnitude  of  the  satellite's  velocity  is  determined  by  the 

i  i       iGM  v 

height  of  its  orbit:    |v|  =  J —  ,  where  G  is  the  Gravitational  Constant  in     /kg  ,  M  is  the 

mass  of  the  earth,  r  is  the  radius  (height)  of  the  orbit,  and  R  is  the  radius  of  the  earth. 
Since  the  satellite's  movement  is  deterministic,  like  the  aircraft  model,  satellite  position 
and  velocity  are  updated  only  when  polled.  Otherwise,  the  satellite  continues  across  the 
sky  and  returns  to  its  starting  point  when  it  reaches  the  ending  footprint  point. 


Figure  5-4.  Satellite  Model  Event  Graph. 


e.   Fixed  Model 

The  simplest  of  all  models  is  the  fixed  sensor.  The  fixed  model's  position 
never  changes  and  its  velocity  is  set  to  zero.  A  call  to  update  a  fixed  object's  position 
literally  results  in  nothing  happening.  Simply,  the  fixed  model  sends  its  information  on 
whenever  polled  and  requires  no  discrete  event  mechanization. 


3. 


Solver 


By  using  a  square  root  information  filter  (SRIF)  in  a  linearized  least-squares 
formulation,  JBSTAFSim's  solver  assumes  that  measurement  errors  are  statistically 
independent  and  distributed  Normal  with  mean  zero  and  input  covariance.  With  some 
modification,  a  weighting  matrix  that  assumes  some  correlation  could  be  implemented, 
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however,    since   little   is   known   about   the   true    distribution    of  TDOA   and    FDOA 
measurement  errors,  an  estimate  of  the  correlation  seems  even  more  unlikely. 

In  order  to  have  an  initial  estimate  of  the  target's  position,  the  linearized  least 
squares  geo-location  solution  is  started  with  the  results  provided  by  a  direct  solution. 
JBSTAFSim  solves  this  set  of  equations  as  described  in  Chapter  IV  using  a  Singular  Value 
Decomposition  routine  adapted  from  Reference  21.  As  discussed  in  Chapter  IV,  after  a 
fix  is  computed,  the  solver  then  uses  the  previous  solution  as  the  starting  position  for  the 
subsequent  solutions  until  a  request  is  made  for  another  direct  solution  to  update  the 
current  solution.  The  periodicity  of  these  updates  may  be  set  by  the  user.  In  an  actual 
Geo-location  System,  direct  solution  updates  might  be  requested  as  a  function  of  many 
complex  variables  and  parameters,  for  instance:  quality  of  the  previous  fix,  the 
dimensionless  shock,  time  since  the  last  fix,  by  request  of  the  system  user,  or  quality  of  the 
received  signal.  While  having  a  set  doctrine  for  computing  direct  solutions  leaves  little 
room  for  qualitative  decision  making  algorithms  that  would  exist  in  the  ARL:UT  system,  it 
is  a  simple  abstraction  for  the  simulation  given  the  assumptions  already  made.  In  addition, 
by  updating  the  target's  initial  "guess"  position,  the  direct  solution  prevents  the  model 
from  losing  track  of  the  target  if  the  previous  linearized  least  squares  solution  was  too 
great  in  error.  The  difficulty  with  any  direct  solution  lies  in  the  dual  solutions  provided  as 
a  result  of  squaring  the  TOA  equations.  As  described  earlier  in  Chapter  IV,  this  ambiguity 
is  generally  resolved  with  some  knowledge  of  where  the  target  is  in  relation  to  the 
observers  or  by  simply  comparing  the  two  solutions  and  throwing  out  the  one  that  is 
obviously  too  distant  from  the  observers.  JBSTAFSim  assumes  that  the  operational 
system  will  have  several  algorithms  for  resolving  this  ambiguity,  and,  thus,  always  picks 
the  correct  starting  point. 


4.  Measurement  Errors 

Related  to  the  solver  error  assumptions,  is  the  actual  error  mechanization  itself. 
JBSTAFSim  measures  precise  TOA's  and  FOA's  in  meters  and  meters  per  second.  These 
values  are  differenced  from  those  of  the  platform  designated  as  the  master.  Thus,  these 
differences  are  the  true  TDOA's  and  FDOA's  that  any  sensor  would  measure  in  the 
absence  of  measurement  noise.  In  order  to  simulate  measurement  noise,  the  truth  values 
are  then  corrupted  by  Normal,  zero  mean  random  values  generated  from  a  distribution 
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defined  by  the  user.  The  user  is  free  to  define  the  standard  deviation  of  errors  across  each 
TDOA/FDOA  measurement.  In  addition,  the  user  is  allowed  to  correlate  the  errors 
between  TDOA/FDOA  measurements. 

The  production  of  measurement  errors  begins  by  first  producing  a  covariance 
matrix  using  the  data  provided  by  the  user.  While  the  user  is  only  concerned  with  the 
correlation  of  measurement  errors  between  those  differences  which  include  the  master  as  a 
sensor  and  not  the  correlation  between  slaves,  for  the  three  dimensional  case,  this 
correlation  must  nevertheless  be  supplied  in  order  to  build  a  covariance  matrix  from  which 
to  generate  random  errors.  Further,  JBSTAFSim  checks  this  value  against  the  partial 
correlation  with  those  errors  associated  with  the  master  to  ensure  that  the  given 
correlation  vector  generates  a  positive-definite  symmetric  covariance  matrix  from  which  to 
generate  errors.  The  covariance  matrix  is  then  factored  into  a  lower  triangular  matrix 
using  the  Cholesky  decomposition  routine.  Next,  JBSTAFSim  uses  SIMKIT's  uniform 
[0,1]  generator  to  produce  Normal  (0,1)  variables  using  the  Box  and  Muller  method 
Calling  the  Normal  (0, 1 )  variables  Z  ,  the  correlated  errors  Xj ,  and  the  lower  triangular 

matrix  C ,  the  measurement  errors  are  produced  from  the  following  algorithm  (Scheuer 
and  Stoller's  method)  from  Reference  23: 

For/  =  l...«,  X,=5XZ.  5-1 

7  =  1 

Equation  5-1  could  be  changed  to  produce  non-zero  mean  random  numbers  by 
simply  adding  the  mean,  fix ,  to  the  sum  on  the  right  hand  side.    Equation  5-1  also  makes 

clear  the  need  for  a  positive  definite  covariance  matrix  since  any  other  matrix  could  not  be 
factored  into  the  matrix  C . 

D.  MODEL  DATA  PROCESSING 

As  discussed  above,  JBSTAFSim  assumes  a  GPS-assisted  geo-location  system  that 
is  functioning  nominally  with  good  communications  between  all  sensors.  JBSTAFSim 
sensors  are  polled  and  fixes  are  computed  every  10  seconds  of  simulation  time.  While  this 
poll  frequency  is  not  set  permanently,  it  was  chosen  to  correspond  in  value  to  that  in  some 
initial    tests    of  the    ARL:UT    system.       JBSTAFSim    never    misses    a    fix    or    loses 


communication  between  sensors  since  these  are  not  characteristics  that  JBSTAFSim  is 
endeavoring  to  simulate. 

The  polling,  fix,  and  result  comparison  process  are  depicted  in  the  flow-chart  in 
Figure  5-5  below. 


Figure  5-5    JBSTAFSim  Flow  Chart 

In  addition  to  scheduling  the  next  poll  sensors  event,  the  poll  sensors  event  causes 
the  target  and  all  sensors  to  update  their  positions  and  velocities.  Sensors,  who  maintain  a 
reference  to  the  target,  compute  their  TOA  (distance  to  target  measured  in  meters)  and 
FOA  (line  of  sight  velocity  measured  in  meters  per  second)  to  the  target.  The  slave 
measurements  are  then  differenced  from  the  master  to  form  N-l  TDOA's/FDOA's.  This 
"truth"  data  is  then  corrupted  by  inducing  measurement  errors  that  are  modeled  by  a 
Normal  distribution  with  zero  mean  and  standard  deviation  and  cross-correlation  input  by 
the  user.  This  data  and  the  sensor  positions  are  then  passed  to  the  solver.  The  solver  uses 
the  sensor  positions  and  the  current  guess  of  the  target's  position  to  form  the  partials 
matrix  for  the  square  root  information  filter  (SRIF).  The  SRIF  then  performs  the  required 
iterations  to  solve  for  the  target's  position.  These  results  are  sent  to  a  data  collector  that 
compares  the  solution  to  the  truth  position  The  data  collector  keeps  running  statistics — 
mean,  variance,  95%  confidence  interval  half-width,  largest  and  smallest  extremes — of  the 
miss  distance,  geometric  dilution  of  precision  (GDOP),  and  size  of  the  area  of  uncertainty . 
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Finally,  the  solver  need  not  be  a  SRTP.  Due  to  JBSTAFSim's  object  oriented 
design,  the  data  returned  by  all  sensors  to  the  master  is  itself  an  object  which  could  be 
processed  by  any  solver.  This  data  object  consists  of  a  collection  of  objects  which  hold 
each  sensor's  position,  velocity,  TO  A  data,  FOA  data,  and  has  methods  to  evaluate  the 
numerical  quantities  which  make  up  the  partials  matrix.  The  master's  data  is  stored  in  the 
first  location  of  this  data  object  vector.  In  order  to  form  the  "DOA's"  from  the  data,  the 
remaining  data  objects  difference  their  information  from  the  master's.  Of  course,  in  a 
linearized  least  squares  formulation,  the  partials  matrix  must  be  evaluated  at  some  point, 
this  point  is  supplied  by  an  object  which  contains  the  current  target  estimate  position  and 
velocity.  When  the  solver  finishes  its  iterations,  it  returns  its  solution  as  this  type  of  an 
object.  Thus,  any  solver  coupled  to  JBSTAFSim  need  only  accept  the  basic  objects  in 
JBSTAFSim  which  define  position  and  velocity  for  any  platform. 

E.  STATISTICAL  RESULT  ESTIMATION  AND  CALCULTAIONS 

Since  the  focus  of  analysis  is  centered  on  the  miss  distance  of  the  geo-location  fix, 
JBSTAFSim  uses  this  parameter's  precision  for  the  terminating  condition.  Miss  distance 
is  simply  the  Euclidean  distance  between  two  points — the  truth  position  and  the  solution. 
Since  both  of  these  points  are  represented  as  objects  in  Java,  the  calculation  of  the  miss 
distance  is  merely  a  matter  of  asking  one  object  the  distance  to  the  other. 

Thus,  JBSTAFSim  seeks  a  precise  estimate  for  the  mean  of  the  miss  distance. 
However,  since  JBSTAFSim  runs  for  no  fixed  sample  size,  the  model  has  no  control  over 
the  confidence-interval  half-length  which  directly  impacts  the  precision  of  the  estimate. 
Therefore,  in  order  to  produce  a  precise  estimate  for  the  mean  miss  distance  and  provide  a 
terminating  condition,  JBSTAFSim  employs  relative  precision.  Using  relative  precision  as 
a  terminating  condition  for  the  estimation  of  the  mean  of  a  random  variable  is  discussed  in 
Appendix  C  and  in  detail  in  Reference  23.  The  advantage  of  using  relative  precision  lies  in 
its  ability  to  provide  a  precise  estimate  of  the  mean  miss  distance  when  the  sample  size  is 
beyond  the  analyst's  control. 

In  addition  to  fix  miss  distance,  JBSTAFSim  keeps  statistical  information  on  the 
geometric  dilution  of  precision  (GDOP)  and  uncertainty  ellipsoid  area  or  volume. 
Ellipsoid  volume  is  calculated  using  the  eigenspace  as  described  in  Chapter  IV,  Equations 
4-19  and  4-20.  JBSTAFSim  employs  a  Jacobi  Transformation  algorithm  which  consists  of 
a  sequence  of  orthogonal  similarity  transformations.     Each  transformation  is  a  plane 


rotation  designed  to  annihilate  one  of  the  off-diagonal  matrix  elements.  While  each 
transformation  may  undo  previously  zeroed  elements,  each  element  continues  to  get 
smaller  and  smaller  until  the  matrix  is  diagonal  to  machine  precision.  [Ref.  21]  With  the 
eigenvalues  of  the  estimate  covariance  matrix,  the  area  or  volume  of  the  ellipsoid  of 
uncertainty  can  be  calculated  as  described  in  Chapter  IV.  GDOP  information  is 
accumulated  to  rate  the  quality  and  consistency  of  the  geometry  of  the  simulation.  GDOP 
calculation  is  a  straightforward  application  of  Equation  4-23. 

These  statistical  results  and  accumulated  data  are  provided  to  the  user  through  the 
model  output  selections  provided  by  the  user  at  the  simulation's  initialization. 

F.  MODEL  INPUT  AND  OUTPUT 

1.  Input 

For  any  simulation  run,  JBSTAFSim  requires  a  large  amount  of  initialization  data. 
The  solver  needs  to  know  how  many  dimensions  and  how  many  sensors  with  which  it  will 
be  dealing.  In  addition  the  user  must  set  the  weights  for  the  TDOA  and  FDOA  data.  In 
practice,  these  weights  would  be  computed  on-line  by  the  system  using  a  hybrid  of  both 
equations  and  algorithms  similar  to  those  in  Chapter  3  and  from  experience  of  the  system's 
performance.  The  solver  can  also  be  set  to  "batch"  mode  in  which  it  combines  the 
information  of  the  previous  fixes  with  the  current  data  to  produce  a  fix.  In  addition  to 
solver  information,  the  user  must  set  the  initialization  data  for  each  sensor  and  target — up 
to  five  total  platforms.  As  discussed  above,  by  inputting  the  error  standard  deviation  and 
error  cross-correlation  between  sensor  measurements,  the  user  must  also  create  a 
distribution  for  the  measurement  errors. 

While  this  amount  of  input  is  extensive,  once  it  is  entered,  the  model  takes  over 
and  requires  no  further  guidance  by  the  user.  In  order  to  make  the  model  initialization 
process  easier,  JBSTAFSim  allows  the  user  to  save  the  input  data  to  a  file  specified  by  the 
user.  With  this  file  of  initialization  data,  the  JBSTAFSim  simulation  run  may  be 
immediately  started  by  simply  selecting  the  input  file.  With  this  construct,  any 
initialization  need  only  be  entered  once;  further  the  data  stored  in  the  input  file  is  relatively 
easy  to  edit  and  could  be  changed  in  any  text  editor. 
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2.  Output 

As  the  user  chooses,  JBSTAFSim  outputs  simulation  results  to  the  screen  and/or  a 
data  file.  If  output  to  a  file  is  selected,  JBSTAFSim  returns  the  results  of  its  simulation  to 
a  file  designated  by  the  user.  This  file  contains  the  statistical  results  of  three  data 
accumulators  for  each  solution  type.  As  described  above,  the  statistical  accumulators 
keep  track  of  the  fix  miss  distance,  GDOP,  and  size  of  the  uncertainty  ellipsoid.  At 
simulation  termination,  the  number  of  data  points,  mean,  variance,  standard  deviation, 
largest  and  smallest  extremes  of  the  data  for  each  accumulator  are  written  to  the  output 
file. 

If  a  graphical  display  has  been  requested  by  the  user,  JBSTAFSim  animates  the 
simulated  movement  of  sensors  and  the  target  while  displaying  the  TDOA  and 
TDOA/FDOA  fixes.  In  addition  a  small  summary  at  the  bottom  of  the  animation  provides 
the  current  fix  miss  distance,  an  updated  average  miss  distance,  and  GDOP  of  the  current 
fix.  As  previously  mentioned,  the  graphical  display  is  in  the  latitude,  longitude,  height 
coordinate  system — although  height  is  not  displayed  in  the  animation.  The  latitude  scale 
on  the  left  of  the  screen  provides  a  quick  reference  for  the  scale  of  the  display.  Each  tenth 
of  a  degree  of  latitude  represents  six  nautical  miles  or  12,000  yards.  Thus,  given  that  the 
coordinates  of  each  sensor  in  the  simulation  is  represented  by  WGS84  Ellipsoid 
coordinates  on  the  order  of  106  -  107  meters  which  translates  to  order  101  latitude  and 
longitude  degrees,  and  is  then  translated  into  integer  screen  pixels  displayed  on  a  small 
scale  (large  area)  chart,  it  should  come  as  no  surprise  that  most  platforms  in  the  display  do 
not  tend  to  move  with  great  speed,  especially  when  fixes  are  being  taken  every  10  seconds 
of  simulated  time. 

Despite  these  facts  and  the  fact  that  the  screen  images  provide  little  analytical 
information,  the  graphical  display  does  provide  a  very  effective  and  quick  means  to 
ascertain  what  is  actually  happening  in  the  model  for  both  the  analyst  and  for  any  audience 
of  the  analysis.  Further,  the  ability  to  provide  such  an  animated  stochastic  simulation 
across  a  variety  of  computer  platforms  must  be  considered  a  true  asset  for  analyst  who 
must  present  the  analysis  to  an  (any)  audience  Yet,  Java  is  able  to  deliver  not  only  on  this 
ability,  but  also  allows  the  analyst  and  audience  to  be  separated  by  thousands  of  miles  with 
only  the  Internet  and  a  Web  browser  between  the  two.  Accordingly,  the  ability  for 
JBSTAFSim  to  function  as  an  applet  is  described  in  the  next  chapter. 
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VI.  JBSTAFSIM  AND  THE  WORLD  WIDE  WEB 

As  its  name  implies,  JBSTAFSim  is  written  in  Java — specifically,  the  first  release 
of  the  Java  Development  Kit  (JDK  1.0.2).  Java  is  an  object-oriented  programming 
language  which  is  currently  gaining  increased  popularity  in  the  programming  world.  One 
of  the  reasons  Java  has  become  so  popular  is  its  run  time  library  which  gives  Java 
bytecode  platform  independence.  The  same  code  (without  recompiling)  can  be  used  on 
Windows  95,  Solaris,  UNIX,  Macintosh,  and  so  on.  Given  the  exponential  growth  of  the 
Internet,  this  is  a  true  advantage,  if  not  necessity,  for  any  programming  language.  As  a 
result  of  this  portability,  JBSTAFSim  can  not  only  run  on  any  of  the  machines  listed  above 
as  an  application,  but  also  can  run  as  an  applet  on  any  machine  using  a  Java  enabled  Web 
browser.  The  JBSTAFSim  applet  requires  no  special  programming,  code,  recompilation, 
or  platform.  It  only  requires  a  Web  browser,  nothing  else.  Therefore,  due  to  Java's 
portability,  JBSTAFSim  can  be  run  and  its  results  examined  instantly  from  anywhere  and 
on  any  machine  with  little  overhead.  Thus,  Java  is  a  full  featured  programming  language 
that  provides  the  ability  to  be  run  from  the  Internet  on  any  machine  via  a  Web  browser. 
There  are  many  more  advantages  to  programming  in  Java,  but  simply  stated,  it  is  a  simple, 
robust,  portable,  high  performance,  and  dynamic  language  which  allows  for  efficient, 
object-oriented,  distributed  code.  For  a  detailed  explanation  of  Java's  design  and 
accomplishments,  Sun  has  published  a  "White  Paper"  which  can  be  found  at 
'  'http://java.  sun.  com/whitePaper/java~whitepaper- 1 .  html. 

Java  programs  that  work  within  a  Web  browser  are  called  applets.  The  idea 
behind  an  applet  is  simple:  "users  download  bytecodes  from  the  Internet  and  run  them  on 
their  own  machines"  [Ref  22].  Some  might  ask  why  use  a  Java  applet  on  a  Web  page  to 
run  programs?  Since  Java  is  a  full  programming  language,  it  has  all  the  power  of  any 
''true"  programming  language  with  the  ability  to  run  dynamically  on  the  host  computer  via 
the  Internet.  Without  Java,  if  you  were  to  design  a  "dynamic"  Web  page  that  could 
respond  to  a  user  or  a  stochastic  simulation,  each  change  in  the  Web  page  would  require 
that  data  be  sent  to  a  CGI  script  on  the  server.  The  script  would  need  to  process  the  data 
and  send  the  results  back  to  the  browser — probably  requiring  the  creation  of  a  new  Web 
page  on  the  fly.  This  is  a  horror  to  program  and  without  a  doubt  will  be  slow.  With  Java, 
once  the  applet's  bytecodes  are  downloaded,  the  power  of  the  programming  language  is 
available  to  run  without  updates  from  the  server.  In  order  to  protect  the  host  computer 
system  from  malicious  bytecodes,  Java  implements  a  "sandbox"  that  restricts  the  applet 
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from  corrupting  memory  outside  the  applet's  process  space.  In  addition,  secure  Web 
browsers,  like  Netscape  Navigator,  or  browsers  that  have  their  security  features  enabled 
by  the  user  (Microsoft  Internet  Explorer)  will  prevent  applets  from  writing  or  reading 
local  files.  Applets  can  include  graphical  user  interfaces,  can  catch  mouse  movements  and 
clicks,  and  can  interpret  text  inputs.  Further,  all  processing  is  performed  by  the  user's 
system;  thus,  the  originating  server  is  not  continually  bombarded  with  hits  for  information 
and  number  crunching.  In  addition,  the  user  is  not  hampered  by  bandwidth  when  running 
an  applet  (although,  bandwidth  could  be  an  issue  downloading  the  applet). [Ref  22] 

The  power  of  Java  and  the  applet  concept  above  can  be  realized  within  the  nature 
and  abilities  of  the  JBSTAFSim  applet  itself.  Quite  simply,  the  only  difference  between  the 
JBSTAFSim  applet  and  JBSTAFSim  application  is  one  Java  class  which  declares  that  it  is 
an  applet  and  will  oversee  the  events  which  occur  in  the  JBSTAFSim  simulation.  With  no 
revisions  or  recompilations,  every  bit  of  code  that  the  simulation  uses  as  an  application  is 
used  by  the  applet.  Thus,  the  full  power  and  stochastic  nature  of  the  JBSTAFSim 
application  can  be  delivered  anywhere  to  a  Web  browser  via  the  Internet  with  precisely  the 
same  code.  An  illustration  of  the  JBSTAFSim  applet  running  in  Netscape  Navigator  is 
shown  in  Figure  6- 1 . 

As  an  example  of  the  power,  flexibility  and  utility  of  this  design,  consider  how  the 
JBSTAFSim  applet  is  currently  configured  to  run  from  the  World  Wide  Web.  A  user, 
somewhere  on  some  operating  system  with  a  Java  enabled  Web  browser  starts  the 
JBSTAFSim  applet  by  visiting  the  appropriate  Web  page.  This  Web  page  currently  exists 
on  a  server  which  is  running  on  a  SPARC  clone  with  NEXT  as  its  operating  system.  The 
JBSTAFSim  bytecodes  are  retrieved  by  the  server  from  another  SPARC  clone  which  uses 
SOLARIS  as  its  operating  system.  The  JBSTAFSim  bytecodes  consist  of  bytecodes 
compiled  on  a  Windows-95  system.  This  code  implements  SIMKIT  bytecodes  which 
were  compiled  on  a  Macintosh.  Further,  the  SIMKIT  bytecodes  implement  bytecodes 
from  the  Java  Generic  Library  that  were  compiled  on  some  system  unknown  to  the  author 
of  JBSTAFSim.  The  server,  ignorant  to  all  of  this,  delivers  all  these  bytecodes  to  the 
user's  system  for  execution  in  the  user's  Web  browser.  The  point  is:  it  is  irrelevant  from 
what  system  the  code  originated,  from  where  the  code  is  coming,  and  to  where  the  code  is 
going  for  execution.  The  bytecodes  need  only  be  delivered  to  the  user's  machine  which 
will  locally  process  the  bytecodes'  information  quickly  and  efficiently — all  bytecodes  run 
the  same  everywhere 
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Figure  6-1   JBSTAFSim  Applet  in  Netscape  Navigator 

For  both  the  analyst  and  the  audience,  this  is  a  powerful  tool.  The  analyst's 
sponsors  can  quickly  inspect  a  product  or  proposal  without  special  software,  hardware,  or 
large  investments  of  time  and  money.  Further,  the  tool  is  extendible  and  can  be  coupled 
with  other  models  regardless  of  the  operating  system  or  environment  in  which  they  were 
originally  created.    In  the  case  of  JBSTAFSim,  in  addition  to  ARLUT  and  faculty  at 
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Naval  Postgraduate  School,  the  applet's  audience  quickly  included  the  Naval  Information 
Warfare  Activity  and  other  Commander,  Naval  Security  Group  offices  and  activities. 
Each  activity  only  requiring  a  Web  browser  and  a  connection  to  the  Internet.  Certainly,  in 
this  light,  JBSTAFSim  could  represent  one  small  component  in  what  has  the  potential  to 
be  a  robust,  high-performance,  and  distributed  modeling  architecture  for  the  DOD  and  its 
community. 
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VH.  SIMULATION  ANALYSIS 

A.  BACKGROUND 

Several  simulation  scenarios  were  conducted  to  evaluate  and  compare  the 
performance  of  both  the  TDOA  and  mixed  TDOA/FDOA  solutions.  Both  two 
dimensional  and  three  dimensional  fix  scenarios  were  conducted.  Within  these  two  cases, 
several  parameters  were  varied  to  quantify  their  effect  on  both  solution  types,  including: 
measurement  error  variance,  the  correlation  between  measurement  errors,  geometry  of 
the  sensor  network,  and  the  frequency  that  the  SRIF  requests  direct  solution  updates. 

B.  TWO  DIMENSIONAL  GEO-LOCATION  SCENARIO  RESULTS 

The  two  dimensional  case  consisted  of  a  sensor  network  of  two  ships  and  a 
satellite  geo-locating  a  ship,  as  depicted  in  Figure  7-1.  This  scenario  was  run  seven  times 
each  for  four  measurement  error  correlation  values.  The  seven  runs  varied:  the 
frequency — continuous,  every  10  fixes,  and  every  20  fixes — of  direct  solution  updates;  the 
TDOA  measurement  error  variance;  and  the  FDOA  measurement  error  variance.  The 
measurement  correlation  values  were:  0.0,  0.1,  0.5,  and  0.9  error  correlation  between  all 
measurements.  The  TDOA  measurement  error  variance  was  examined  for  the  values  of 
602  and  302  meters.  The  FDOA  measurement  error  variance  was  examined  for  the  values 
of  0.22  meters  per  second  and  0.12  meters  per  second.  Sensor  to  target  geometry  was 
mediocre  with  GDOP  averaging  1.5  with  minimum  and  maximum  values  of  1.3  and  2.6, 
respectively.  For  reference,  the  values  of  602  meters  and  0.22  meters  per  second  for  the 
respective  TDOA  and  FDOA  measurement  error  variances  are  approximately  equivalent 
to  those  empirically  experienced  in  the  ARL:UT  system. 

1.  Geo-Location  Accuracy 

Results  from  these  simulation  runs  are  depicted  in  Figure  7-2  through  Figure  7-17 
with  some  summarized  information  provided  Table  7-1.  The  distance  from  the  solution's 
position  to  the  target's  actual  position — the  miss  distance — is  the  primary  measure  of 
effectiveness  (MOE),  while  the  area  of  the  2-sigma  uncertainty  ellipse  is  a  secondary 
MOE.  Thus,  the  primary  MOE  is  a  measure  of  the  solution's  accuracy  while  the 
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Figure  7-1  Two  Dimensional  Scenario  Sensor  Allocation 


secondary  MOE  is  a  measure  of  the  solution's  ability  to  rate  its  precision.   In  terms  of  the 
miss  distance,  the  primary  MOE,  Figures  7-2  through  7-8  show  that  the  mixed 
TDOA/FDOA  solution  maintains  first  order  stochastic  dominance  for  each  simulation  run 
which  varied  the  fix  update  frequency  in  addition  to  the  cases  where  the  TDOA  and 
FDOA  measurement  error  variances  were  contrasted. 

a.  Fix  Updates 

In  order  to  prevent  JB ST AF Sim's  solution  from  diverging  due  to  a  poor 
target  estimate,  direct  solutions  are  provided  at  a  regular  update  interval.    This  interval 


46 


was  varied  from  continuous  updates,  to  updates  after  every  ten  fixes,  and  for  updates  after 
every  20  fixes.  Clearly,  as  seen  in  Figure  7-2,  the  continuous  direct  solution  updates  are 
degrading  the  accuracy  of  both  solutions,  but  particularly  so  for  the  TDOA  solution. 
These  continuous  update  distributions  are  so  poor,  they  are  dismissed  from  further 
analysis.  Notice  that  while  the  TDOA  solution  and  mixed  T/FDOA  solution  are  more 
similar  in  the  "10  fix  update"  and  "20  fix  update"  cases,  the  "10  fix  update"  scenario 
offers  the  best  shape  to  this  distribution.  This  can  be  seen  in  the  relative  location  of  the 
point  in  the  distribution  where  the  curve  makes  its  "right  hand  turn"  in  the  upper  quantiles. 
Notice  that  this  location  tends  to  be  higher  and  more  to  the  left  in  the  "10  fix  update" 
scenario. 


Comparison  of  TDOA  and  FDOA  Miss  Distance  Quantiles 
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Figure  7-2  Comparison  of  TDOA  and  T/FDOA  Miss  Distance  Quantiles 
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Comparison  of  TDOA  and  FDOA  Miss  Distance  Quantiles 
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Figure  7-3  Comparison  of  TDOA  and  T/FDOA  Miss  Distance  Quantiles 
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b.  Solution  Type  Performance 

In  each  scenario  above,  the  mixed  T/FDOA  case  dominates  the  TDOA 
solution.  The  T/FDOA  solution's  dominance  is  largest  in  the  upper  portion  of  the 
distribution  after  the  characteristic  curve  as  noted  above.  In  order  to  show  these  contrasts 
more  clearly,  Figure  7-9,  below,  plots  some  of  the  scenario  results  together  In  this  figure, 
a  horizontal  line  representing  the  0.95th  quantile  is  drawn  to  intersect  the  solution 
distribution  curves.  This  line  intersects  the  mixed  T/FDOA  solution,  with  measurement 
error  variances  of  302  m  and  0.22  m/s,  100  meters  to  the  left  of  the  TDOA  solution  with  a 
measurement  error  variance  of  302  m.  Thus,  the  is  a  significant  difference  between  the 
two  solution  types  in  the  upper  areas  of  their  distributions. 


Comparison  of  TDOA  and  T/FDOA  Miss  Distance  Quantiles 
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Figure  7-9  Comparison  of  Several  TDOA  and  T/FDOA  Miss  Distance  Quantiles 

Even  more  interesting,  notice  that  the  mixed  solution  curve  identified 
above  (error  variances  of  302  m  and  0.22  m/s)  outperforms  a  mixed  solution  with  an 
improved  FDOA  error  variance  of  0.1 2  m/s  and  a  degraded  TDOA  measurement  error 
variance  of  602  m.  Further,  note  that  at  this  quantile,  this  second  mixed  solution  curve 
(error  variances  of  60"  m  and  0.12  m/s)  has  started  to  dominate  the  TDOA  solution  with 
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an  improved  TDOA  measurement  error  variance  of  302  m.  Thus,  this  T/FDOA  solution 
outperforms  that  of  the  TDOA  solution  in  the  upper  quantiles  where  either  measurement 
error  or  geometry  is  degrading  the  accuracy  of  the  fix.  Therefore,  despite  a  larger  error 
variance  in  the  TDOA  measurement  errors,  a  mixed  T/FDOA  solution  provides  more 
accurate  fixes  when  facing  situations  that  are  more  difficult  for  the  TDOA  solution  to 
solve  by  itself.  Thus,  the  designers  or  users  of  a  system  such  as  this  would  desire  a  larger 
reduction  the  TDOA  measurement  error  variance  as  compared  to  that  in  the  FDOA 
measurement  even  though  it  is  the  FDOA  information  which  will  improve  the  accuracy  of 
the  fix. 

c.    Upper  Quantile  Performance 

The  most  striking  feature  of  these  curves  is  the  sharp,  near  right-angle  turn 
that  the  distributions  take  at  the  upper  quantiles.  The  effect  of  rotating  any  of  the  quantile 
plots  90  degrees  to  the  left  so  that  the  actual  quantiles  are  on  the  X-ax\s  make  the  "right- 
angle  turn"  feature  of  the  distribution  appear  even  more  pronounced.  With  this 
understanding  of  the  shape  of  the  distribution,  gaining  an  intuition  of  the  difference 
between  the  two  solutions  that  transcends  first  order  stochastic  dominance  can  be  gleaned. 
As  pointed  out,  the  area  in  which  the  two  solution  distributions  differ  the  most  can  be  seen 
in  the  separation  between  the  curves  in  the  higher  quantiles.  In  this  light,  the  mixed 
T/FDOA  solution's  dominance  can  be  applied  to  the  outliers  of  the  distribution.  In  the 
cases  where  noise  or  geometry  degrade  the  TDOA  solution,  the  additional  information 
from  the  FDOA  measurements  increases  the  accuracy  of  the  T/FDOA  solution  over  that 
of  TDOA  solution.  Thus,  it  is  especially  in  these  cases  of  bad  geometry  or  large  noise  that 
the  mixed  solution  outperforms  the  TDOA  only  solution.  Unfortunately,  at  each 
individual  data  point,  it  is  unclear  from  these  results  whether  it  is  the  geometry  or  the 
measurement  error  which  causes  the  size  of  the  error  in  the  location  of  the  target  Two 
interesting  plots  of  the  mixed  T/FDOA  performance  distribution  are  shown  in  Figures 
7-10  and  7-11  below.  Notice  in  Figure  7-10  that  the  lower  FDOA  measurement  error 
variance  leads  to  better  results  as  would  be  expected.  In  addition,  holding  the  FDOA 
measurement  error  variance  at  0.2"  m/s  and  decreasing  the  TDOA  measurement  error 
variance  leads  to  a  distribution  which  dominates  the  other  two  distributions.  These  are  not 
unexpected  results,  however,  Figure  7-11  turns  these  results  around.  The  distribution 
with  the  worst  measurement  error  variance  dominates    Notice  that  Figure  7-10  displays 
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T/FDOA  Miss  Distance  Quantiles  by  Sigma  TDOA  and  FDOA 
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Figure  7-10  Comparison  of  TDOA  and  T/FDOA  Miss  Distance  Quantiles 
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distributions  where  the  measurement  error  correlation  is  high — 0.9;  yet,  Figure  7-11 
portrays  the  results  from  0.0  measurement  error  correlation.  It  is  difficult  to  describe  why 
measurements  with  higher  error  correlation  should  outperform  those  with  no  correlation. 
Further,  it  is  not  obvious  that  measurements  with  a  larger  error  distribution  should 
outperform  those  with  better  measurement  data.  Yet,  this  is  exactly  what  these  figures 
indicate.  These  results  lead  to  two  avenues  of  discussion — geometry  and  uncertainty. 

If  the  results  of  Figure  7-10  and  7-1 1  were  functions  of  geometry,  it  is  not 
clear  from  the  GDOP  measurements  the  degree  to  which  geometry  is  degrading  the 
accuracy  of  the  fix.  Recall  from  Chapter  IV  that  GDOP  is  not  really  a  measure  of 
geometry  but  a  description  of  the  diagonal  elements  of  the  estimate  covariance  and 
weighting  matrices.  It  is  more  a  measure  related  to  what  the  system  thinks  might  happen 
versus  that  of  what  is  actually  happening  in  the  system.  In  Reference  13,  Chaffee  and 
Abel,  the  authors  of  JBSTAFSim's  direct  solution  theory,  discuss  how  the  null  space  of 
the  direct  solution  relates  to  the  two  roots  of  the  direct  solution.  Perhaps,  the  null  space 
could  provide  some  information  with  regard  to  the  geometry  of  the  observers  and  their 
target,  however  this  type  of  analysis  might  only  describe  the  geometry  of  the  TDOA 
surfaces.  It  is  more  than  conceivable  that  the  ability  of  the  mixed  T/FDOA  solution  to  be 
more  accurate  in  the  highest  quantiles  of  the  miss  distance  distribution — the  outliers — can 
be  attributed  to  cases  where  the  geometry  of  the  TDOA  surfaces  is  poor  yet  the  geometry 
of  the  FDO  A  surfaces  is  sufficient  enough  to  improve  the  mixed  solution.  In  view  of  the 
shape  of  the  miss  distance  distributions  in  the  upper  quantiles,  this  theory  could  describe 
the  sharp  "right  hand  turn"  as  described  above  The  solution  to  rating  the  geometry  of  a 
network  of  observers  in  relation  to  their  target  is  beyond  the  scope  of  this  thesis;  however, 
JBSTAFSim  would  be  the  perfect  tool  to  empirically  experiment  with  this  theme  of 
geometry.   The  second  issue,  fix  uncertainty,  deserves  a  section  unto  itself  to  describe. 

2.  Uncertainty  of  G^eo-locations 

If  geometry  is  not  all  or  only  part  of  the  answer  to  the  unexplained  improvement  in 
the  accuracy  of  the  solution  in  the  presence  of  correlation,  then,  perhaps,  the  uncertainty 
of  the  estimate  may  provide  another  avenue  to  help  explore  this  issue.  Figures  7-12 
through  7-17  display  the  distribution  of  the  area  of  the  two-sigma  uncertainty  ellipse 
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associated  with  both  solution  types  for  the  continuous,  10  fix,  and  20  fix  update  scenarios. 
In  these  graphs,  the  standard  span  (the  area  inside  the  "whiskers")  of  the  distribution  is 
drawn  at  1.5  times  the  Inter-Quartile  Range  beyond  the  quartiles.  All  points  outside  the 
standard  span  are  drawn  individually.  With  this  construction,  it  is  easy  to  view  the  density 
of  the  outliers. 

When  the  TDOA  solution  is  compared  with  that  of  the  mixed  T/FDOA  solution, 
the  main  difference  between  the  two  plots  is  the  density  of  the  outliers.  Another  way  of 
indicating  the  larger  spread  of  the  T/FDOA  uncertainty  ellipse  size  is  through  the 
comparison  of  the  standard  deviation  of  this  area.  These  quantities  are  summarized  in 
Table  7-1  below. 


Error 
Correla- 

Continuous 

Update 

10  Fix  Update 

20  Fix  Update 
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Sigma 

Sigma 

Sigma 

tion. 

TDOA 

T/FDOA 

TDOA 

T/FDOA 

TDOA 

T/FDOA 
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13197.1 
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10962 
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13175.6 

55483.4 
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4901 

9659 
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55616.6 

5643.1 

14718.5 

4703 

9658 

0.9 

13708.7 

55613.3 

5576.7 

14228.6 

4610 

11007 

Table  7-1  Distribution  of  TDOA  and  T/FDOA  Ellipse  Size  Standard  Deviations 

First,  notice  that  the  ellipse  size  standard  deviations  for  the  mixed  T/FDOA 
solutions  are  much  larger  than  those  of  the  TDOA  solutions  in  all  cases.  This  is  confirmed 
by  the  outlier  spread  and  densities  in  the  graphs  below.  Second,  notice  that  in  the  10  Fix 
Update  scenario  (where  the  accuracy  of  the  0.0  correlation  was  compared  to  that  of  the 
0.9  correlation  in  Figures  7-10  and  7-11  above)  the  solution  which  produced  more 
accurate  results — 0.9  correlation — has  a  larger  standard  deviation  of  ellipse  area. 
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Figure  7-13  Comparison  of  T/FDOA  Ellipse  Areas  by  Correlation 
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TDOA  2-Sigma  Ellipse  Area  Distributions  by  Meas.  Error  Corr. 
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Figure  7-14  Comparison  of  TDOA  Ellipse  Areas  by  Correlation 
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TDOA  2-Sigma  Ellipse  Area  Distributions  by  Meas.  Error  Corr. 
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Now  it  is  possible  to  tie  geometry,  error  correlation,  accuracy,  and  uncertainty 
together.  Clearly,  the  T/FDOA  solution  dominates  in  accuracy;  however  it  tends  to  have  a 
larger  spread  of  uncertainty  with  regards  to  its  results.  Yet,  comparing  the  median  and 
inter-quartile  ranges  of  the  uncertainty  ellipse  sizes  for  the  two  solutions,  they  are  very 
similar.  Thus,  the  mixed  T/FDOA  solution  is  able  to  transcend  errors  induced  by 
geometry  and  measurement  errors,  possibly  induced  by  correlation,  and,  thus,  provide 
more  accurate  fixes — especially  in  the  upper  quantiles  of  the  accuracy  distributions.  The 
price  for  such  accuracy  is  reflected  in  the  uncertainty  of  the  estimate.  This  should  not  be 
upsetting,  and,  in  fact,  proves  the  robust  nature  of  the  Square  Root  Information  Filter. 
Since  the  SRJF  reports  each  solution  with  no  idea  with  regard  to  its  accuracy  (accuracy 
requires  knowledge  of  the  true  target  location),  it  can  only  rate  each  fix  by  the  solution 
estimate  covariance.  Therefore,  even  though  the  SRTF  provides  a  more  accurate  fix,  it 
correctly  reports  that  measurement  errors,  error  correlation,  or  geometry  are  degrading  its 
confidence  in  the  solution  that  it  has  provided. 

C.  THREE  DIMENSIONAL  GEO-LOCATION  SCENARIO  RESULTS 

The  three  dimensional  scenarios  included  a  base  line  "good"  geometry  scenario,  to 
show  the  effects  of  error  correlation,  and  the  comparison  of  a  "good"  and  "bad" 
geometry  cases  to  show  the  effects  of  sensor  allocation  on  geometry.  No  dominance  of 
either  solution  type  for  the  three-dimensional  case  is  shown. 

Setting  up  scenarios  to  compute  three  dimensional  fixes  is  very  difficult.  This 
difficulty  is  due  mostly  to  the  geometry  of  the  scenario.  A  three  dimensional  fix  requires 
four  observers  in  order  to  form  three  differences  for  the  evenly  determined  case.  The 
addition  of  the  height  dimension,  of  course,  expands  both  the  partials  matrix  (from  a  2x2 
to  a  3x3  matrix  for  the  TDOA  solution  and  from  a  4x4  to  a  6x6  matrix  for  the  mixed 
T/FDOA  solution)  and  the  data  vector  (from  2x1  to  3x1  and  from  4x1  to  6x1  for  TDOA 
and  T/FDOA  formulations,  respectively)  in  the  linearized  least  squares  formulation.  As  a 
result  each  difference  pair  between  slave  and  master  must  allocate  its  position  information 
from  the  partials  matrix  to  its  respective  TDOA  and/or  FDOA  measurement(s).  The 
difficulty  in  geometry  arises  from  the  formulation  of  partial  matrices  from  these  difference 
pairs  that  are  near  singular.  While  this  is  straightforward  in  theory,  it  is  difficult  to  explain 
concisely  to  a  warrior  allocating  his  or  her  assets  to  perform  a  passive  geo-location 
Further,  the  possibility  of  the  construction  of  a  near  singular  matrix  is  not  as  readily 
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apparent  when  viewing  the  relative  three-dimensional  geometry  as  compared  to  that  of  the 
two-dimensional  geometry.  Finally,  and  most  profoundly  for  an  operational  system,  this 
problem  is  further  exacerbated  when  one  considers  that  in  practice,  the  warrior  may  have 
little  information  with  regard  to  the  target's  true  position  and,  thus,  defining  the  geometry 
a  priori  becomes  an  even  more  difficult  task.  Yet,  this  does  not  mean  three-dimensional 
fixes  need  be  abandoned.  The  scenarios  which  depict  "good  geometry"  provide  excellent 
fixes — many  with  median  miss  distances  well  below  100  meters — that  outperform  poor 
geometry  two-dimensional  fixes.  The  main  difficulty  with  a  three-dimensional  fix,  then, 
may  be  providing  the  warrior  with  the  means  to  optimally  allocate  his  or  her  assets. 

1.  Good  Geometry  Results 

The  good  geometry  scenario  consisted  of  a  sensor  network  of  two  ships,  an 
aircraft,  and  a  fixed  sensor  targeting  a  ship.  A  visual  depiction  of  the  scenario  is  provided 
in  figure  7-18.  This  scenario  was  run  for  20  measurement  error  correlation  values  from 
-0.9  to  0.95.  For  the  purposes  of  identifying  the  correlation  combinations,  a  simulation 
run  labeled  with  a  negative  correlation  value  actually  means  that  within  the  three  TDOA 
measurements  and  within  the  three  FDOA  measurements,  two  measurements  are 
correlated  negatively  by  a  magnitude  equal  to  the  positive  correlation  of  the  other  pair  of 
measurements.  For  example,  if  measurements  1  and  2  were  negatively  correlated,  then 
measurements  1  and  3  would  be  positively  correlated  while  measurements  2  and  3  would 
be  negatively  correlated. 

The  results  of  this  scenario  are  multifaceted.  First,  restricting  the  sensors  and  the 
target  into  very  tight  boxes  ensured  that  the  geometry  of  the  TDOA  measurements  would 
be  almost  always  perfect.  Thus,  the  TDOA  measurements  dominated  the  solution  with 
regard  to  the  position  of  the  target.  Therefore,  the  identical  results  of  both  the  TDOA  and 
T/FDOA  solutions  proves  the  stability  and  accuracy  of  the  mixed  solution's  square  root 
information  filter — a  particularly  gratifying  result  given  the  formulation,  scaling,  and 
un-scaling  of  the  linearized  least  squares  equation  as  described  in  Chapter  IV  and 
Appendix  B. 

The  surprising  result,  as  depicted  in  Figure  7-19,  is  the  improvement  in  the 
accuracy  of  the  solution  as  the  correlation  between  measurement  errors  increases. 
However,  upon  inspecting  Figure  7-20,  the  opposite  is  true  of  the  distribution  of  the 
volume  of  the  uncertainty  ellipses.     Given  the  results  from  the  two-dimensional  case 
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Java-based  Stochastic  TDOA/FDOA  Simulation 
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Figure  7-18  Geometry  of  Basic  Scenario 

discussed  above,  this  is  no  longer  a  surprise.  While  the  SRjT  has  no  idea  that  it  is 
producing  more  accurate  results  in  the  highly  correlated  noise  environment,  it  is  robust 
enough  to  realize  that  the  measurement  errors  which  it  is  encountering  are  deviating  from 
what  the  weighting  matrix  had  told  it  to  expect  The  SRJF  produces  the  best  fix  it  can; 
yet,  it  honestly  reports,  in  the  uncertainty  of  the  fix,  that  the  measurement  errors  are  not 
indicative  of  what  it  had  expected 
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TDOA  and  T/FDOA  Miss  Distance  Distributions  by  Meas.  Error  Correlatioi 
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Figure  7- 1 9  Comparison  of  Miss  Distance  Distributions  by  Correlation 


TDOA  and  T/FDOA  2-Sigma  Ellipse  Volume  Distns.  by  Meas.  Error  Corr. 
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Figure  7-20  Comparison  of  Ellipse  Volumes  by  Correlation 
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Comparison  of  TDOA  and  T/FDOA  Miss  Distance  Quantiles  by  Error  Corr 
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Figure  7-21  Comparison  of  Miss  Distance  Distributions  by  Correlation 
2.  Comparison  of  Good  and  Bad  Geometry  Scenarios 

As  described  in  the  introduction  to  this  section,  the  geometry  of  the  observers 
relative  to  their  target  is  crucial  for  an  accurate  three-dimensional  solution.  Figures  7-22 
and  7-23  show  the  configurations  of  the  "bad  geometry"  and  "good  geometry"  scenarios 
for  this  illustration,  respectively.  Figures  7-24  and  7-25  show  the  quantile  distributions  of 
these  two  scenarios. 

Examining  the  distribution  of  sensor  platforms  in  Figure  7-22,  it  is  not  clear — even 
with  the  position  of  the  target  known — that  this  allocation  of  sensor  assets  would  lead  to  a 
"bad  geometry"  scenario,  yet,  the  GDOP  in  this  scenario  averaged  222  with  minimum  and 
maximum  values  of  149  and  1472,  respectively.  Compare  this  to  the  "good  geometry 
scenario  in  Figure  7-23  whose  GDOP  averaged  0.8962  with  minimum  and  maximum 
values  of  0.8874  and  0.91277,  respectively.  Without  this  data  from  the  simulation  runs,  it 
would  have  been  impossible  to  predict  these  results.  Neither  case  appears  exceptionally 
threatening  from  the  geometry  standpoint;  yet,  it  will  be  the  scenario's  geometry  that 
directly  impacts  the  accuracy  of  the  geo-location  solution. 
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-.Java-based  Stochastic  TDuA/FOOA  Simulation 
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Figure  7-22  Bad  Geometry  Case 
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Java-based  Stochastic  TDOA/FDOA  Simulation 
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Figure  7-23  Good  Geometry  Case 
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Good  Geometry  TDOA  and  T/FDOA  Miss  Distance  Quantifies 
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Figure  7-24  Good  Geometry  Miss  Distance  Distribution 
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Figure  7-25  Bad  Geometry  Miss  Distance  Distribution 

Upon  examining  Figures  7-24  and  7-25,  the  possible  impact  of  poor  geometry  on 
the  accuracy  of  the  three-dimensional  solution  becomes  clear     An  explanation  of  what 


66 


makes  the  "good  geometry"  case  "good"  is  required.  In  terms  of  sensors,  the  difference 
between  the  cases  is  the  replacement  of  one  of  the  aircraft  observers  from  the  "bad"  case 
with  a  satellite  in  the  "good"  case.  Given  its  altitude,  the  satellite  is  able  to  create  a  better 
three-dimensional  surface  when  its  position  is  differenced  from  that  of  the  master  As  a 
result,  the  satellite  "sees"  height  better  than  that  of  the  aircraft  who  may  as  well  be  on  the 
surface  of  the  earth  in  comparison  to  the  satellite's  height.  The  partials  matrix,  then, 
contains  more  useful  information  with  regard  to  the  third  dimension.  Thus,  the  surface 
defined  by  the  difference  of  the  satellite's  and  the  master's  information  is  better  suited  to 
complement  the  other  surfaces  formed  by  the  other  observers  and,  therefore,  can  provide  a 
more  accurate  fix  with  this  "good  geometry"  of  measurement  surfaces. 

The  second  obvious  feature  when  examining  Figures  7-24  and  7-25,  is  that,  like 
the  first  good  geometry  case  in  the  previous  section,  both  the  TDOA  and  T/FDOA  miss 
distance  distributions  are  identical.  In  this  analysis,  the  author  was  not  able  to  produce  a 
case  where  the  T/FDOA  solution  stochastically  dominated  the  TDOA  solution.  While  the 
construction  of  such  a  case  may  be  possible,  performance  of  both  solutions  would  most 
likely  be  highly  dominated  by  the  quality  of  sensor  geometry  versus  that  of  sensor 
measurements.  Further,  while  individual  cases  can  be  shown  where  a  single  TDOA  or 
mixed  T/FDOA  solution  was  obviously  better  than  that  of  the  other,  the  author  was  not 
able  to  show  that  this  lead  to  the  stochastic  dominance  of  either  solution  type.  Thus,  the 
solution  to  improving  three-dimensional  fixes  does  not  appear  to  be  in  the  augmentation  of 
FDOA  information,  however,  FDOA  information  did  not  on  the  whole  appear  to  degrade 
the  distribution  of  the  three-dimensional  solution  either.  The  consequences  from  these 
results  lead  to  two  avenues  of  discussion:   solution  formulation  and  sensor  selection. 

A  solution  to  improving  the  three-dimensional  fix  may  be  simply  collapsing  to  two- 
dimensions.  This  analysis  has  limited  itself  to  the  exploration  of  solutions  in  the  evenly 
determined  case.  Yet,  it  is  likely  that  most  overdetermined  two-dimensional  solutions  will 
dominate  or  at  least  do  as  well  as  those  in  the  evenly  determined  case.  In  the 
overdetermined  case,  it  has  not  been  shown  to  what  extent  FDOA  information  could 
improve  or  degrade  the  solution.  Likewise,  the  quality  of  a  three-dimensional  fix  from  an 
overdetermined  sensor  network  with  more  than  four  observers  has  not  been  shown.  Yet, 
having  a  plethora  of  sensors  to  distribute  at  will  seems  overly  optimistic  and  would  not 
suit  the  warrior  who  desires  to  know  a  lower  bound  on  the  performance  of  his  or  her 
equipment.  Yet,  these  are  both  important  analysises  which  should  be  performed  to 
understand  the  operational  capabilities  of  the  true  system. 

67 


On  the  other  hand,  if  the  utility  for  a  three-dimensional  fix  is  high,  then  the  warrior 
who  needs  such  high  quality  information  should  be  concerned  with  the  selection  of  robust 
three-dimensional  assets.  In  particular,  in  this  analysis  it  has  been  shown  that  the  satellite 
is  an  excellent  sensor  to  assist  in  observing  height.  Yet,  it  was  also  seen  in  two- 
dimensional  simulation  runs  that  in  the  rare  case  when  a  satellite  passes  directly  overhead 
the  target,  it  can  say  little  about  the  target's  two  dimensional  location  since  the  satellite 
only  "sees"  height.  Again,  the  proper  allocation  of  assets  sometimes  may  be  out  of  the 
warrior's  control.  Given  that  satellite  coverage  may  not  always  be  available,  another 
three-dimensionally  robust- sensor  might  be  one  of  the  prototype  unmanned  aerial  vehicles 
(UAV's)  that  is  designed  to  operate  at  very  high  altitudes  for  long  periods  of  time — for 
example,  the  Tier  11+  Conventional  High  Altitude  Endurance  (HAE)  UAV  (Global  Hawk) 
or  Tier  Ill-Low  Observable  HAE  UAV  (Dark  Star)  both  of  which  advertise  surveillance 
altitudes  greater  than  18  kilometers  and  12  kilometers,  respectively.  Note  that  the  first 
three-dimensional  "good  geometry"  scenario  only  used  one  aircraft.  Further,  both  UAV's 
and  satellites  may  be  less  obtrusive  ways  of  improving  the  sensor  to  target  geometry  when 
airspace  restrictions  or  the  requirement  for  the  concealment  of  sensors  dictate  that  the 
traditional  sensors  would  not  suffice. 

Like  most  decisions  a  military  commander  must  make,  the  selection  and  allocation 
of  sensors  to  perform  passive  geo-location  is  not  a  straightforward  choice.  Yet,  further 
analysises  of  like  those  suggested  above  could  assist  the  warrior  in  this  decision  process. 
Again,  a  good  tool  for  exploring  any  of  these  issues  would  be  a  slightly  modified  version 
ofJBSTAFSim. 
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Vffl.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.  CONCLUSIONS 

A  simulation  based  methodology  was  used  to  explore  the  differences  between  geo- 
location  solutions  computed  using  both  time-difference-of-arrival  (TDOA)  and  mixed 
TDOA  and  frequency-difference-of-arrival  (FDOA)  techniques  using  a  Multi-platform, 
GPS-assisted  architecture.  This  methodology  assumed  no  one  method  to  measure  either 
the  TDOA  or  FDOA;  yet,  allowed  for  the  stochastic  modeling  of  the  errors  in  such 
measurements  in  such  a  general  way  that  they  could  apply  to  any  TDOA  or  FDOA 
measurement  process.  Likewise,  the  platforms  which  measure  such  quantities  were 
modeled  in  abstract  ways,  yet  retained  each  platform's  most  basic  behavior.  In  addition, 
the  simulation  portrayed  joint  platforms  operating  together  in  a  joint  operating 
environment.  To  enhance  the  numerical  veracity  of  the  solutions,  the  simulation 
implemented  a  numerically  stable  Square  Root  Information  Filter  (SRIF)  which  received 
initial  target  estimates  based  on  a  numerically  stable  direct  solution  implementing  the 
Singular  Value  Decomposition  routine.  These  numerically  robust  methods  were  used  to 
provide  the  most  numerically  stable  and  accurate  solutions  possible  in  the  face  of  poor 
sensor  geometry  and  computer  round-off  error.  Data  taken  during  the  simulation 
included:  the  solution  miss  distance  from  truth,  the  area  or  volume  of  the  2-sigma  ellipse 
centered  around  the  solution,  and  the  geometric  dilution  of  precision  (GDOP).  The 
analysis  focused  on  the  accuracy  of  the  solutions  represented  by  the  fix  miss  distance,  and 
the  solution  uncertainty,  represented  by  the  volume  or  area  of  the  2-sigma  ellipse  around 
the  fix. 

1.  Two  Dimensional  Conclusions 

In  the  two-dimensional  problem,  it  was  shown  that  a  mixed  T/FDOA  solution 
maintains  first  order  stochastic  dominance  over  the  TDOA  solution  in  mediocre  geometry. 
The  price  the  mixed  solution  pays  for  this  increased  accuracy  is  a  large  increase  in  the 
spread  of  its  fix  uncertainty;  however,  both  solutions  maintained  roughly  the  same 
partitioning  of  fix  uncertainty  around  the  center  of  their  respective  uncertainty 
distributions.  Therefore,  increasing  the  spread  of  the  uncertainty  distribution  remains  a 
small  price  to  pay  for  the  increased  accuracy  of  the  solution.     The  largest  difference 
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between  fix  accuracy  was  seen  in  the  upper  quantiles  of  miss  distances  where  the  mixed 
T/FDOA  solution  significantly  improved  the  solution  fix  accuracy  over  that  of  the  TDOA 
solution.  Thus,  the  FDOA  surfaces  appear  to  substantially  improve  the  solution  accuracy, 
especially  when  the  TDOA  surfaces  would  lead  to  poor  fixes. 

In  addition  to  the  contrasts  between  solution  types,  it  was  also  shown  that  while 
the  direct  solution  prevents  the  linearized  least  squares  solution  from  diverging,  the 
amount  to  which  the  SRIF  was  updated  by  the  direct  solution  significantly  impacted  the 
accuracy  of  both  solution  types.  In  particular,  solutions  which  used  the  direct  solution  as 
an  initial  estimate  for  every  fix  were  so  poor,  their  results  were  disregarded. 

2.  Three  Dimensional  Conclusions 

As  an  exploration  into  the  impact  of  correlated  errors,  measurement  error 
correlation  was  varied  across  20  values  in  a  good  geometry  scenario.  While  both 
solutions  produced  more  accurate  results  in  the  higher,  positively  correlated  cases,  the 
solutions  also  indicated  more  uncertainty  in  these  fixes.  Again,  this  confirmed  the  robust 
nature  of  the  SRIF  formulation.  These  results  are  useful  since  the  actual  system  would 
implement  sensors  taking  measurements  in  a  common  operating  environment  with 
common  noise  characteristics,  thus,  most  likely  creating  similar  measurement  errors. 

The  three-dimensional  simulations  did  not  show  stochastic  dominance  of  either 
solution  type.  While  some  simulations  showed  differences  in  the  accuracy  of  individual 
solution  type  fixes,  neither  solution  type  was  shown  to  stochastically  dominate  the  other. 
Results  which  produced  identical  accuracy  distributions  for  both  the  TDOA  and  mixed 
T/FDOA  solutions  across  several  scenarios  proved  the  stability  of  the  algorithm  which 
solves  the  mixed  T/FDOA  solution.  The  inference  drawn  from  these  simulations  was  that 
geometry,  more  than  solution  and  measurement  type,  impacts  the  accuracy  of  the  three- 
dimensional  fix  solution.  Therefore,  the  choice  of  sensors  which  can  more  optimally 
measure  three-dimensional  information  should  be  the  first  concern  for  those  who  require 
such  fixes.  In  particular,  it  was  shown  that  sensors  with  high  altitude  characteristics  are 
well-suited  for  this  task. 
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B.  RECOMMENDATIONS 

1.  System  Use  and  Design 

The  Multi-Platform,  GPS-Assisted  TDOA  Geo-location  System  as  designed  by  the 
Applied  Research  Laboratories,  University  of  Texas,  Austin  is  a  robust  and  potentially 
highly  functional  passive  geo-location  system  for  both  the  Navy  and  the  Department  of 
Defense.  FDOA  information  can  be  safely  added  to  the  geo-location  formulation  to 
augment  the  TDOA  information  for  the  solution  of  the  target's  location.  In  the  two- 
dimensional  solution  case,  the  FDOA  information  can  significantly  improve  the  accuracy 
of  the  geo-location.  In  the  three  dimensional  case,  the  additional  information  does  not 
appear  to  degrade  the  system's  performance  in  any  significant  manner. 

The  linearized  least  squares  formulation  in  the  SRIF  appears  to  be  a  very  stable 
and  robust  formulation  for  solving  the  geo-location  problem.  The  system  which 
implements  this  formulation  should  implement  some  direct  solution  or  other  means  to 
provide  an  roughly  accurate  estimate  of  the  target's  location  to  start  the  SRIF's  solution 
iterations.  While  the  actual  system  would  not  be  computing  the  large  number  of  fixes  for 
one  target  that  were  required  in  this  stochastic  analysis,  a  means  to  periodically  update  the 
SRIF  with  a  new  target  estimate  should  be  implemented.  Further,  the  actual  system  would 
need  to  request  this  information  based  on  signal  characteristics,  computation  of  the 
dimensionless  shock,  or  some  combination  of  both  solver  and  signal  properties.  In 
addition,  such  an  algorithm  might  be  able  to  re-weight  the  partials  matrix  in  the  presence 
of  correlated  noise.  The  ability  to  chose  a  more  optimal  measurement  covariance  matrix 
would  significantly  improve  the  fix  estimate  covariance  and,  thus,  decrease  the  uncertainty 
of  the  estimate  of  the  target's  location. 

An  avenue  to  improving  three-dimensional  fixes  was  shown  to  be  through  the  use 
of  platforms  which  can  measure  height  information  more  reliably.  This  consequence  leads 
to  two  recommendations.  First,  the  ARL:UT  system  should  be  incorporated  and  tested  on 
satellite  systems  and  as  a  potential  payload  for  prototype  UAV  platforms.  Second,  since 
this  matter  is  at  its  crux  a  geometry  issue,  a  new  method  to  rate  the  geometry  of  the  sensor 
network,  other  than  GDOP,  should  be  pursued.  For  the  warrior,  in  particular,  such  a 
method  should  be  able  to  assist  in  rating  the  sensor  network's  geometry  given  possible 
locations  of  the  target.   Such  a  method  would  allow  the  commander  to  allocate  his  or  her 
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assets  in  an  pseudo-optimal  way  given  either  possible  locations  of  the  target  or  given  no  a 
priori  information  with  regard  to  the  target's  location. 

2.  Java-Based  Stochastic  TDOA  and  FDOA  Simulation 

JBSTAFSim  proved  itself  as  a  robust  stochastic  simulation  for  the  analyses 
presented  in  Chapter  VII.  In  order  to  study  any  of  the  proposed  areas  above,  modified 
versions  of  JBSTAFSim  would  provide  the  perfect  tool  for  the  analysis  of  any  of  these 
issues.  In  addition,  JBSTAFSim  could  be  modified  so  that  specific  platforms  or  platform 
characteristics,  for  example  a  particular  type  of  satellite  orbit,  could  be  evaluated  for  its 
ability  to  contribute  to  the  geo-location  solution.  In  addition,  due  to  time  and  scope 
limitations,  JBSTAFSim' s  "batch  mode"  option  for  its  solver  was  never  stochastically 
tested  with  a  static  target.  JBSTAFSim  requires  more  programming  so  that  it  can  make 
qualitative  decisions  with  regard  to  the  frequency  with  which  it  resets  its  batch  mode 
solver. 

Investigating  solution  types  in  both  the  overdetermined  two-dimensional  and  three- 
dimensional  cases  should  be  undertaken.  While  it  is  likely  that  these  solutions  would 
dominate  the  evenly  determined  ones,  it  is  not  clear  to  what  extent  FDOA  information 
might  affect  these  solutions.  This  would  be  a  particularly  interesting  study  given  that 
overdetermined,  "bad  geometry"  three-dimensional  fixes  could  potentially  be  collapsed  to 
more  accurate  overdetermined  two-dimensional  fixes.  With  some  effort  JBSTAFSim 
could  be  modified  to  analyze  this  issues. 

Next,  JBSTAFSim  has  modeled  both  the  TDOA  and  FDOA  measurement  errors 
as  Normally  distributed  random  variables.  The  Normal  distribution  was  implemented  for 
its  ease  in  the  ability  to  generate  correlated  random  variables  on  a  computer  and,  as 
presented  in  Chapter  HI,  due  to  the  fact  that  any  measurement  error  distribution  is  in  itself 
a  function  of  many  random  variables  which  change  over  a  variety  of  system  configurations 
and  environmental  characteristics.  Quite  simply,  the  true  distribution  is  neither  known  nor 
can  it  be  known;  yet  it  can  be  modeled.  To  that  end,  modeling  the  measurement  errors 
with  another  distribution  would  prove  an  interesting  test  for  the  SRIF  and  would  be 
relatively  straightforward  to  implement  in  JBSTAFSim.  Investigating  the  solution 
accuracy  and  uncertainty  using  either  distributions  that  are  either  "heavy  tailed"  or  skewed 
would  be  worthy  endeavors.  Likewise,  extending  JBSTAFSim's  solver  to  create  a  non- 
diagonal  measurement  error  covariance  matrix  might  prove  useful  in  any  of  these  analyses. 
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Of  the  many  possibilities  that  exist  for  the  extension  of  JBSTAFSim,  one  final 
possibility  could  be  of  particular  interest.  Given  its  design  and  network  aware  capabilities, 
JBSTAFSim's  existing  modular  design  could  be  modified  so  that  it  receives  actual  TDOA 
and  FDOA  measurements  and  sensor  information  via  a  network.  JBSTAFSim's  solver 
would  then  compute  the  geo-location  and  display  the  results  on  a  local  computer  or  back 
to  the  network  to  be  displayed  on  another  host's  computer  using  only  a  Web  browser. 
While  achieving  this  vision  requires  the  resolution  of  many  security  issues,  Public  Key 
Encryption  and  secure  sockets  capabilities  within  the  Java  language  itself  are  already  being 
implemented  in  the  latest  version  of  Java.  Further,  this  concept  again  shows  the  "flat" 
architecture  that  is  innate  in  models  created  in  Java.  With  very  little  cost  to  the  user  or 
developer,  rich  capabilities  are  immediately  available  to  the  model  which  can  function 
across  many  platforms  in  a  very  robust  nature. 
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APPENDIX  A.  SEQUENTIAL  SQUARE  ROOT  DATA  PROCESSING 

ALGORITHM  AND  THEORY 


The  purpose  of  this  appendix  is  to  describe  both  the  underlying  theory  behind  the 
square  root  information  filter  (SRIF)  and  the  sequential  data  processing  algorithm  which 
exploits  this  theory.  A  basic  understanding  of  solving  linear  least  squares  problems  is 
assumed,  but  is  not  essential.  The  primary  reference  of  this  section  is  Reference  16  thus, 
notation  will  be  presented  as  similarly  as  possible.  In  essence  this  section  will  provide  a 
quick  summary  reference  of  the  SRIF  algorithm  as  described  in  detail  in  Reference  16. 
For  the  most  part,  all  applicable  ideas  and  examples  presented  here  will  be  identical  to 
those  in  Reference  16. 

A.  CONTEXT  OF  THE  PROBLEM  AND  REVIEW  OF  LEAST  SQUARES  [REF. 
16, 14-16] 

With  regard  to  geolocation  in  a  F/TDOA  problem,  the  most  basic  sequential  linear 
least  squares  system  of  equations  is  defined  as: 


cF 

not 


X  —  JCo)  =z~F{xo) 


A-Al 


Xo 


Where  ^/L-  represents  the  TDOA  partials  or  mixed  F/TDOA  partials  evaluated 
at  the  state  update  x0 .   This  MxN  matrix  is  multiplied  by  the  state  step  N  vector  [x-  x0) 

which  is  approximately  equal  to  the  M  data  vector  z  minus  the  TDOA  or  F/TDOA 
functional  M  vector  evaluated  at  x0 . 
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From  Chapter  IV  and  for  reference,  the  TDOA/FDOA  equations  are  provided 


below. 


Observation  equations: 


c-TDOA1J(r)  =  Rl')-Ru) 


f;lFDOAv(r,v)  = 


Av£>      Av£> 


Measurement  covariance  matrix: 


M  = 


c2a2 

*-    u  TDOA 


0 

fTC<7 


2„2  _2 

FDOA 


(a0)-a0)) 


Av 


,tO) 


Av 


«(o 


^0) 


/? 


(0 


Linearized  least  squares  equation: 


(a0,-fl(i)) 


r  c-TDOAl}  \ 
—  FDOAv 


,Av<<>  -  Av^ 


*■„ 


Direction  cosines  are  denoted  by  a ,  and  i?  is  the  observer-target  separation.  In 
the  state  vector,  dx  and  dv  refer  to  the  offset  from  the  position  and  velocity,  respectively, 
of  the  nominal  state,  and  the  subscripts  p  and  _L  refer  to  components  parallel  and 
perpendicular,  respectively,  to  the  line  of  sight.  The  symbols  c  and  fT  refer  to  the  speed 
of  light  and  the  frequency  of  the  transmitter.  Of  course,  the  dimension  of  the  state  here  is 
6,  not  2  as  it  appears,  because  all  the  elements  on  the  left  side  of  the  equation  are 
vectors.  [Ref  19] 

For  simplicity,  let  us  now  consider  A-Al  as  the  set  of  equations: 


Ax  +  v  =  z 


A-A2 
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Here  the  partials  matrix  in  A-Al  is  represented  by  the  MxN  matrix  A ;  the  state  step  by  the 
vector  x ;  and  the  functional  and  data  difference  by  is  represented  by  the  vector  z  .  We  also 
introduce  v  to  represent  the  M  vector  of  observation  errors.  The  least  squares  solution  to 
minimize  the  mean  square  observation  error  is  defined  as: 

•/  =  X>0)2  =  ^  A-A3 

Or,  in  terms  of  x  : 

J(x)  =  (z-Ax)T(z-Ax)  A-A4 

Thus,  J(x)  is  non-negative  and  quadratic  in  the  components  of  x  so  that  the  necessary  and 
sufficient  conditions  for  a  minimum  require  that  x  satisfy  the  normal  equations: 

ATAx=ATz  A-A5 

A-A5  leads  to  the  well-known  least-squares  solution  xh : 

Xls=(ATAyATz  A-A6 

Statistically,  if  we  assume,  for  now,  that  the  observation  errors,  v ,  are  random 
variables  that  have  been  normalized  so  that: 

E(v)  =  0  and  E(vTv)  =  IM  A-A7 

A-A5  then  becomes: 

AT  Axh  =  ATz  =  ATAx  +  ATv  A-A8 

A-A8  allows  us  to  make  two  important  conclusions  about  our  least  squares 
solution.   First,  that  E(xfa)  =  x ,  thus  it  xb  is  an  unbiased  estimate.    Second,  the  estimate 

error  covariance  can  be  evaluated  by: 
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(ata)e[(Xis  -x\Xls  -x)T](ATA)  =  ATE(vvT)A  A-A9 

A-A9  simplifies  so  that  we  get  a  recognizable  estimate  for  P^,  defined  as  estimate 
error  covariance: 

P*  =  ^{xu-xfa-xf]  =  (ATA)-'  A-A10 


Where  P .   =  A  A  is  defined  as  the  information  matrix,  A  . 


B.  A  PRIORI  INFORMATION  [REF.  16, 16-17] 

Now  suppose  that  in  addition  to  the  linear  system  A-A2,  there  exists  an  a  priori 
unbiased  estimate  of  x  and  A  ,  defined  as  x  and  A  .  This  information  can  be  included  in 
the  least  squares  solution  by  changing  A-A4  to: 

J,(x)  =  (x -  x)  A(x -x)+(z-  Ax)  (z  -  Ax)  A-Bl 

Again,  like  A-A5,  the  minimizing  argument  for  «/,(*)  ,  now  define  as  xh ,  is  defined 
by  satisfying  the  normal  equations: 

(X  +  ATA)xh  =  7uc  +  ATz  A-B2 

From  A-B2  and  being  reminded  that  (x  -  x)  and  v  both  have  mean  zero,  a  similar 
construct  to  A-A8  can  be  written  to  prove  that  x^is  unbiased.   Further,  assuming  that  the 

observation  "noise"  and  the  a  priori  estimate  errors  are  independent  (£'[(x-x)v:r]  =  Oj, 
then  like  A-A9  and  A-A10  we  obtain  the  estimate  error  covariance: 

pa=4(*-4x*-%)r]=(A+^r  a-b3 

Notice  how  A-B2  and  A-B3  reduce  to  A-A5  and  A-A10  when  A  =  0,  i.e.  there  is 
no  a  priori  information. 
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C.  A  PRIORI  INFORMATION  AS  ADDITIONAL  OBSERVATIONS  [REF.  16, 
17-18) 

Given  that  factorization  of  a  positive  semi-definite  matrix  is  always  possible,  factor 
A  so  that  A  =  RTR  .  RT  is  called  the  square  root  of  A  .  Square  root  matrices  are  not 
unique,  but  are  related  by  orthogonal  transformations.  For  example,  if  5,  and  S^are 
square  roots  of  a  positive  semi-definite  matrix,  then  P  =  SxSf  =  S2Sl ,  there  exists  a 
matrix  J  such  that  S2  =  ^Tand  TTT  -  TTT  =  /  . 

With  this  factorization  in  mind,  let  z  -  Rx  and  write: 

(x  -  x)T  A(jc  -  x)  =  (z  -  Rx)T(z  -  Rx)  A-C 1 

Compare  these  results  to  J(x)  in  A-A4  and  one  can  glean  that  a  priori  estimate- 

covariance  information  can  be  interpreted  as  additional  observations  in  the  least  squares 
solution.    Therefore,  J,(x)  in  A-Bl  can  be  viewed  as  simply  applying  the  least  squares 

functional  J{x)  to  the  augmented  system: 


A-C2 


A-C2  is  a  recursive  result  since  the  previous  estimate  and  covariance  are  combined 
with  the  new  data  to  form  an  updated  estimate  and  covariance.  Numerically  significant  is 
that  the  estimate  covariance  results  do  not  depend  on  which  square  root  was  used — 
allowing  the  selection  of  square  root  matrices  that  are  simple  to  compute  and  enhance 
numerical  accuracy. 

D.  HOUSEHOLDER  ORTHOGONAL  TRANSFORMATIONS  [REF.  16,  57-62] 

In  section  C,  we  made  a  brief  reference  to  orthogonal  transformations.  SRIF 
depends  on  the  Householder  orthogonal  transformation.  An  understanding  of  this 
transformation  is  vital  in  order  to  understand  the  SRIF  mechanization. 

First,  a  review  of  orthogonal  matrices.  A  matrix  T  is  orthogonal  if  TT  =  I . 
Four  properties  useful  to  know  about  orthogonal  transformations  are: 
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z 

~R~ 

V 

z 

— 

A 

x  + 

V 

A.  If  Tx  and  T2  are  orthogonal,  so  is  T}T2 . 

b.  r1  =  tt  and  rrr  =  /. 

C.  For  any  M  vector  j/,  ||7y||  =  \y\  where  ||-||   is  the  Euclidean  norm. 

D.  If  vis  an  M  vector  of  random  variables  with  E(v)  =  0  and  E[vvT)  =  I , 
then  v  -  Tv  has  the  same  properties  as  v ,  namely:  E{v)  =  0  and  E\vvT)=  I . 

Recall  from  A-A4  that  for  MxN A ,  J{x)  =  {z-  Ax)T(z  -  Ax)  =  \\z  -  Axf .    Define 
T  to  be  a  MxM  orthogonal  matrix  and,  using  Property  C  from  above,  we  have: 


J(x)  =  \\T(Z-Axf  =\\(Tz)-(TA)x\\ 


A-Dl 


Since,    xfa    is   independent   of    T,    T    can   be   chosen   so   that    (TA)    has   a 
computationally  attractive  form  such  as: 


TA  = 


m-n 


A-D2 


Where  R  is  an  upper  triangular  matrix.  If  Tz  is  partitioned, 


Tz  = 


I" 
\m-n 


A-D3 


then  J(x)  can  be  written  as: 


J(x)  = 


T,       II2  II  I'2 

RX\\     +\\Zr 


A-D4 


From  A-D4,  one  can  see  that  the  minimizing  xu  must  satisfy  Rx  =  z,   and  that 

/       \        (I      112 

J\xls)  =  \\z2\\  .    These  results  have  been  proved  to  be  less  susceptible  to  errors  due  to 

computer  round-off  in  References  17  and  18.  From  above,  we  see  the  applicability  of 
orthogonal  transformations  to  the  least  squares  problem.  Let  us  now  focus  on  an 
orthogonal  transformation  that  will  assist  us  in  the  results  from  section  C. 
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Householder  transformations  are  matrix  representations  corresponding  to  the 
geometric  notion  of  a  reflection.  Let  a  bea  nonzero  vector  and  let  U±  be  the  plane 
perpendicular  to  «.  If  v  is  an  arbitrary  vector, 


y  =  (yTu)u  +  v 


A-D5 


then  the  reflection  of  y  in  the  plane  U±  is  equal  to: 


yr  =  -(yTu)u  +  v 


A-D6 


where  u  is  a  unit  vector  in  the  direction  of  u  and  v  is  that  part  of  y  that  is  orthogonal  to 
u  .   See  Figure  A-Cl.  Eliminating  v  from  both  equations  gives: 


yr  =  y  -  2(yTu)u  =  (l-  puuT)y 


A-D7 


where  /?  =  z~lT  ■  The  matrix  T  =1  -  /3uuT  is  the  elementary  Householder  transformation. 

Ml 

Figure  A-Cl  shows  the  geometry  of  the  Householder  transformation. 


Figure  A-Cl 

Five  properties  of  this  Householder  transformation,  Tu ,  are  important  to  us. 
1-     Tu  =  TuT ,  Tu  is  symmetric. 

2.     T^  =  / ,    Tu   is  idempotent.     (Combined  with  Property  1   shows  that    Tu   is 
orthogonal.) 


81 


3.  If  u\j)  =  0  then  [Tj^j)  =  y(j) ,  rather,  if  the/th  component  of  u  is  zero,  then 
Tu  leaves  theyth  component  of  y  unchanged. 

4.  If  u±y  then  Tuy  =  y. 

5  Tuy  =  y-yu,  where  y  =  2yTu/uTu .  This  is  practically  a  time  and  storage 
saver  since  forming  and  storing  Tu  prior  to  computing  Tuy  requires  more 
computation  than  the  direct  evaluation  using  this  property.  If  M  is  large  this 
could  be  particularly  helpful. 

6.  Let  a  =  sgn(^(l)) •  (yTy)  2  and  define  u(\)  =  y{\)  +  a ,  u(j)  =  y(j),j  >  1 . 
Then  Tuy  =  -ae}  and  2/uTu  =  l/(aw(l)) . 

Property  6  allows  us  to  choose  a  direction  of  u  so  that  yr  =  Tuy  lies  along  el 
which  gives  7^y  =  (a,0,...,0)  .  This  is  the  first  step  in  matrix  triangularization.  Since 
Property  C  in  our  orthogonal  transformation  review  assures  us  that  the  length  of  y  is 

(t   \n 
y  y)    .   The  direction  of  u    can  be  obtained  from 

Property  5  above:     u  =  constyy  -  ae}J ,  which  explains  Property  6  again- uTu  =  2ou(\)  . 

Refer  back  to  Figure  A-Cl;  note  yr  -y  is  orthogonal  to  Ul  and  so  is  parallel  to  u  .  Also 

note,  in  Property  6,  that  the  sign  of  o  is  chosen  to  maximize  [«(l)|  this  assists  in  reducing 

numerical  problems  when  Tu  is  applied  to  the  columns  of  A  as  \/cru(\) . 

If  Property  6  is  applied  to  an  MxN  matrix  A  (with  y  as  its  first  column)  the  first 

step  of  a  matrix  triangularization  results. 


T,A=  ) 

m-  \\ 


n-\ 


A 


A-D8 


Note  that  via  Property  6,  s  and  A  are  computed  directly  from  A  and  the  matrix 
Tu  is  implicit — never  calculated.  Repeated  application  of  Property  6  and  the  results  of 

A-C8  to  the  columns  of  A  results  in  a  series  of  orthogonal  transformations  which 
combine  into  one  orthogonal  transformation,  T  (by  Property  A  of  orthogonal 
transformations)  and  can  be  written  as: 
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TA  = 


0 


tf, 


a. 


A-D9 


sk  at 


*k+\ 


E.  THE  SRIF  DATA  PROCESSING  ALGORITHM  [REF.  16,  69-72] 

The  last  four  sections  now  culminate  into  the  SRIF  Data  Processing  Algorithm. 
This  method  is  reputed  to  be  more  accurate  (less  susceptible  to  the  effects  of  computer 
round-off  errors)  and  stable  (accumulated  round-off  errors  will  not  cause  the  algorithm  to 
diverge)  than  other  conventional  non-square  root  methods. 

Begin  by  supposing  that  we  have  an  a  priori  array  from  A-Cl  [Rz\  corresponding 

to  the  equation  z  =  Rx  +  v  where  v  has  zero  mean  and  unity  covariance  and  R  is  NxN. 
We  are  interested  in  constructing  the  least  squares  solution  to  the  a  priori  data  equation 
and  the  new  measurements  z  -  Ax  +  v .  In  section  C  equation  A-C2,  the  least  squares 
solution  was  described  as: 


z 

~R 

V 

z 

— 

A 

x  + 

V 

A-El 


In  section  D,  this  was  seen  to  be  equivalent  to  finding  the  least  squares  solution  to: 


~R~ 

x  =  T 

z 

-T 

V 

A 

z 

V 

A-E2 


Further,  section  D  showed  that  T  can  be  chosen  so  that: 


R 

A 


R 

0 


A-E3 
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where  R  is  an  upper  triangular  matrix  and  we  set: 


z 

= 

z 

and  T 

V 

= 

V 

z 

e 

V 

_Ve_ 

Then  A-E2  becomes: 


Rx  =  z-v 

and 
0  =  e-v„ 


A-E4 
A-E5 


A-E4  is  in  the  form  of  a  data  equation  and  can  now  act  as  the  a  priori  data 
equation  that  will  be  combined  with  the  next  set  of  observations.  The  e  in  A-E5  is  the 
error  in  the  least  squares  fit;  recall  from  A-D4: 


J(x)  =  \\Rx-z\\  +\\e\\ 


A-E6 


Thus,  ||e|    is  the  sum  of  squares  residual  error  corresponding  to  the  least  squares  solution. 
The  algorithm  then  requires  the  construction  of  the  augmented  information 

and  then  can  be  expressed  as  follows  (letting  the  a  priori  information  array 


R 
A 


array 

equal  to  [^z0]) 


3h 

L4 


z 


R. 


J  =  \,...,m 


A-E7 


where    j/^i0 1  =  [^20 1    and    T}    orthogonal  to  triangularize  the  matrix 


R 


7-1 

A. 


Thus, 


Xj  =  R~xz]  is  the  least  squares  estimate  corresponding  to  the  a  priori  information  and  the 
measurement  set,  and  P}  =  R~^R~T  is  the  covariance  of  this  estimate.  Note  that  because 
Rj  is  triangular,  the  calculation  of  these  two  estimates  can  be  easily  accomplished  using 
back  substitution  methods. 
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F.  WHITENING  OBSERVATION  ERRORS  [REF.  16,  47-49] 

Early  in  section  A,  equation  A-A7  asserted  thalE[vTv)  =  IM.  We  now  consider 
the  case  where  E\  vT  v)  =  Pv  where  Pv  is  positive  definite.  We  know,  then,  that  we  can 
write  this  matrix  as: 


P.  =  L..LI 


A-Fl 


where  Lv  is  a  lower  triangular  square  root  of  Pv  .  Now  we  can  multiply  the  observations 
by  L~J  and  the  observation  set  will  have  unit  covariance  observation  error.  Thus  the 
familiar  Ax  +  v  =  z  is  transformed  to  Ax  +  v  =z  with: 


-i 


-i. 


=  L~lz,    A=VvlA,    v  =  Z> 


In  the  TDOA  and  TDOA/FDOA  case  we  assume  that  the  observations  are  uncorrelated, 
but  the  rth  observation  error  has  variance  a,2 .  Thus,  L~J  =  Diag{l  I  a, , . .  .1  /  <jm) . 

Therefore,  our  augmented  matrix  takes  on  the  form: 


R         z 

v}.a  l:1z 


A-Fl 


As  a  side  note,  Reference  16  includes  a  description  of  the  case  where  observations  are 
correlated. 
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APPENDIX  B.  MIXED  TDOA/FDOA  SCALING 

The  purpose  of  this  appendix  is  to  discuss  the  scaling  of  mixed  TDOA  and  FDOA 
partials,  covariance  and  data  in  the  linear  least  squares  formulation.  What  follows  is  a 
summary  of  and  the  examples  directly  from  an  unpublished  work  by  Dr.  Brian  Tolman  at 
ARL:UT,  Reference  19. 

First,  consider  the  usual  linear  least  squares  formulation  with  measurement 
covariance  M  as  described  in  Ref.  19: 

AX  =  B  B_! 

X  =  (ATM~lA)'1  AtM~xB  B-2 

c  =  (atm-'a)-1  b"3 

where, 

A  =  —  B_4 

G*  x0  B-5 

B  =  D-F{X0) 

Dr.  Tolman  then  considers  two  arbitrary  diagonal  matrices  to  re-scale  the  problem. 

S1,  =  diagonal  (mxm)  0.5 

S2  =  diagonal{nxri) 

A  -  o,  AS2  g_g 

X=S~1X  B_9 

B  =  S,B  B-10 

M=SXMSX  B-ll 

The  problem  and  its  solution  then  become 
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AX-B  B-12 

X={ATM~1Ayl  ATM~lB  =  S;i(ATM-1A)~1  ATM~lB  =  S~'X 

C  =  {ATM-X2y  =  S;x{ATM-'AyS-x  =  S^CS^1  B-14 

[Note  that  if  S  is  diagonal,  the  effect  of  multiplying  on  the  left  (SA)  is  to  multiply  row  (i) 
by  S(i,i),  and  the  effect  of  multiplying  on  the  right  is  to  multiply  column  (j)  by  S(j,j).] 

The  scaling  by  S,  can  be  thought  of  as  multiplying  each  observation  equation  by  a 

separate  constant,  as  when,  for  example,  the  units  of  the  TDOA  equation  are  changed 
from  seconds  to  meters.  The  scaling  by  S2  can  be  thought  of  as  changing  the  units  of  the 

elements  of  the  state  vector  (which  of  course  also  changes  the  solution  covariance).  The 
consistency  of  the  above  equations  shows  that  the  linear  problem  may  be  arbitrarily  scaled 
without  modifying  the  solution.    In  non-linear  problems  as  well,  the  scaling  by  S,  clearly 

acts  in  the  same  way,  effectively  multiplying  each  observation  equation  by  a  (perhaps 
different)  constant.    However,  the  scaling  by  S2  cannot  work  for  non-linear  problems, 

except  in  special  cases,  because  there  the  matrix  A  is  in  fact  a  matrix  of  partial  derivatives 

2  JF{X)  _*<*),  B-15 

ax       ex    2 

and  the  second  equality  here  would  hold  only  if  F(X)  were  linear  in  X.  The  TDOA 
equations,  however,  are  a  special  case  (see  below),  and  the  state  vector  can  be  scaled, 
with  one  important  caveat.  [Ref  19] 

Numerical  analysis  indicates  that  the  ideal  situation  is  to  have  a  matrix  A  which  is 
"equilibrated,"  meaning  that 

EKhZKhco/wr-V*  B-16 


That  is,  minimizing  the  numerical  error  incurred  while  solving  the  least  squares  equation 
requires  that  the  sum  of  absolute  values  along  every  row  and  every  column  has  about  the 
same  magnitude.  While  this  is  possible  in  principle,  in  practice  it  is  very  difficult  to  find 
scaling  matrices  which  will  do  this  in  the  general  case.  [Ref.  1 9] 
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Consider  in  detail  the  scaling  of  the  mixed  F/TDOA  problem  as  described  in  Reference  19. 


Observation  equations: 


c-WOAJr)  =  R0)-RU) 


f;1FDOA1J(r,v)  = 


Av« 


Av<>> 


B-17 
B-18 


Measurement  covariance  matrix: 


M  = 


c2<72 

C    U  TDOA 


0 
fT-2c2a 


2 
FDOA 


B-19 


Linearized  least  squares  equation: 


(tt0)-a0)) 
Av[j)     Av{'r 


R 


O) 


R 


(0 


(aU)-a0)) 


\dv) 


(  c-TDOA,,  } 


-FDOA,, 


'  R(>)-RU)  y 
AvJ?  -  Av£>, 


B-20 


Direction  cosines  are  denoted  by  a ,  and  R  is  the  observer-target  separation.  In 
the  state  vector,  dx  and  civ  refer  to  the  offset  from  the  position  and  velocity,  respectively, 
of  the  nominal  state,  and  the  subscripts  p   and   J_   refer  to  components  parallel  and 

perpendicular,  respectively,  to  the  line  of  sight.  Note  the  scaling  of  these  equations.  The 
TDOA  equation  has  units  (meters)  and  typically  these  values  are  of  order  105  or  more. 
The  FDOA  equation  has  units  (m/s)  and  typical  values  are  of  order  1 .  Therefore  the 
relative  scaling  is  of  order  lO"5,  and  may  approach  much  smaller  values.  This  is 
particularly  true  when  the  problem  begins  to  approach  singularity.  Thus  there  is  a 
significant  difference  of  scale  within  the  mixed  TDOA/FDOA  problem.  Also  note  that  the 
TDOA,  FDOA  and  mixed  problems,  although  non-linear,  can  be  scaled  in  the  sense  of 
scaling  the  state  vector,  as  long  as  all  three  components  of  position  and  all  three 
components  of  velocity  are  scaled  alike.  To  see  this  it  is  enough  to  note  that  changing  the 
scale  of  the  three  position  (velocity)  components,  uniformly,  does  not  change  the  physics 
of  the  problem,  it  simply  changes  the  units  of  length  (or  time  or  both).[Ref.  19] 
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Thus  the  F/TDOA  problem  may  be  scaled  in  the  following  way.  Choose  to  scale 
the  problem  with  4  constants,  one  each  for  the  three  position  coordinates  together,  the 
three  velocity  components,  the  TDOA  data  and  the  FDOA  data.  The  scaling  matrices  are 


$  = 


S2  = 


{L'x 

0) 

I  o    jr!J 

(I  <$\ 

U     vj 

B-21 
B-22 


where  each  element  represents  a  diagonal  block  with  constant  elements;  the  blocks  in  S2 
have  dimension  3,  and  the  blocks  in  S,  cover  the  TDOA  data  (L)  and  FDOA  data  (V). 
Then  the  least  squares  equation  becomes 


(aU)-aU))j  0 

Li 

(0\  /  v 


RU)         Rd 


J- TDO^ 
FDO^ 

{    fTVlc    ) 


\r"-rU))/l 


Av 


,-tO) 


V     V 


V 


B-23 


with  measurement  covariance 


M  = 


c2a2 

u    u  TDOA 


L2 

0 


FDOA 


f2V2lc' 


B-24 


Examination  of  these  equations  allows  one  to  identify  an  appropriate  value  for 
each  scaling  constant  so  as  to  make  the  terms  of  order  unity,  as  much  as  is  possible.  The 
length  scale  L  should  be  the  scale  of  the  TDOA  data,  which  is  of  order  the  typical 
observer-observer  separation  distance.  The  other  length  scale,  / ,  is  of  order  R ,  the 
observer-target  separation.  If  these  length  scales  are  equal  the  TDOA  problem  is  well 
conditioned  and  the  partials  matrix  has  all  terms  of  order  unity.  If  they  are  not,  the  TDOA 
problem  is  difficult  to  solve  (the  geometry  is  poor),  but  this  is  a  problem  which  scaling  is 
not  going  to  address  anyway.  In  any  case  it  is  probably  simplest  and  best  just  to  choose  L 
=  I  =  an  average  or  other  compromise  of  the  two  length  scales.  [Ref.  19] 
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In  the  FDOA  portion,  V  apparently  needs  to  be  equal  to  two  velocity  scales,  both 
the  scale  of  the  FDOA  data  (right  hand  side  of  the  equation),  which  has  the  scale  of  the 
typical  observer-observer  relative  velocity  along  the  target  line-of-sight,  and  also  the  scale 
of  the  typical  observer-target  relative  velocity  perpendicular  to  the  line-of-sight  (left  hand 
side).  Again,  if  these  scales  are  equal  the  problem  is  well  conditioned,  and  if  they  are  not, 
there  is  nothing  to  do  about  it.  The  second  velocity  scale  is  arbitrary,  since  it  appears  only 
once,  in  the  partials  matrix;  it  is  probably  best  used  to  cancel  the  scaling  of  the  direction 
cosines,  i.e.  v  =  V .  With  this  scaling  in  place,  clearly  every  equation  is  dimensionless, 
and  every  term  has  a  magnitude  of  the  order  of  unity  (or  at  least  close).  Of  particular 
importance  is  the  partials  matrix,  which  is  nicely  "equilibrated"  This  should  have  a  direct 
effect  on  the  numerical  stability  and  accuracy  of  the  estimation  process.  [Ref.  19] 
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APPENDIX  C.  RELATIVE  PRECISION 

Relative  precision  provides  a  precise  estimate  of  the  mean  when  the  sample  size  of 
the  data  is  unknown.  Relative  precision  requires  that  we  continue  to  sample  our 
simulation  results  until  we  are  certain,  at  some  level  of  confidence,  a ,  that  the  current 
estimate  has  a  given  relative  error,  y  .  Another  way  of  saying  this  would  be  that  we  are 
(l  -  a)%  confident  that  the  percentage  error  in  our  estimate,  the  empirical  mean,  is  lOOy 
percent.  From  Reference  23,  this  probabilistic  relationship  is  defined  as: 


( 


X-m 


M 


\ 


<r 


(l-a) 


C-l 


In  order  to  satisfy  this  relationship,  JBSTAFSim  samples  the  miss  distance  until  there  have 
been  at  least  30  observations  and  the  confidence  interval  half- width  of  the  estimate 
diveded  by  the  empirical  mean  is  less  than  or  equal  the  percentage  error  in  our  estimate. 
Rather,  calling  the  empirical  variance  S2   and  Zk  the  lOOAth  quantile  for  a  standard 

Normal  random  variable,  this  terminating  condition  is  defined  as: 


C-2 


Thus,  the  left  hand  side  of  C-2  represents  an  estimate  of  the  actual  relative  error. 
Once  y  is  set,  the  simulation  continues  until  this  estimate  satisfies  the  above  inequality 

which  means  the  point  estimate  X  has  a  relative  error  of  at  most  y/(\-y)  with 
probability  of  approximately  (l-a).  The  derivation  of  C-l  and  C-2  can  be  found  in 
Reference  23 . 

Finally,  it  should  be  noted  that  our  estimate,  X ,  is  actually  distributed  with  a 
Student-/  distribution  with  n  degrees  of  freedom,  however  since  n  will  be  very  large 
(greater  than  1000),  the  quantile  from  the  standard  Normal  distribution  is  more  than 
adequate  to  calculate  the  estimate  confidence  interval  half-width. 
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